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Abstract 

In  the  last  few  years  meshless  methods  for  numerically  solving  partial 
differential  equations  came  into  the  focus  of  interest,  especially  in  the 
engineering  community.  This  class  of  methods  was  essentially  stimulated 
by  difficulties  related  to  mesh  generation.  Mesh  generation  is  delicate  in 
many  situations,  e.g.,  when  the  domain  has  complicated  geometry;  when 
the  mesh  changes  with  time,  as  in  crack  propagation,  and  remeshing  is 
required  at  each  time  step;  when  a  Lagrangian  formulation  is  employed, 
especially  with  non-linear  PDE’s.  In  addition,  a  need  to  have  flexibility 
in  the  selection  of  approximating  functions  {e.g.,  the  flexibility  to  use 
non-polynomial  approximating  functions),  played  a  significant  role  in  the 
development  of  meshless  methods.  There  are  many  recent  papers,  and  two 
books,  on  meshless  methods;  most  of  them  are  of  engineering  character, 
without  any  mathematical  analysis. 

In  this  paper  we  address  meshless  methods  and  the  closely  related 
generalized  finite  element  methods  for  solving  linear  elliptic  equations, 
using  variational  principles.  We  give  a  unified  mathematical  theory  with 
proofs,  briefly  address  implementational  aspects,  present  illustrative  nu¬ 
merical  examples,  and  provide  a  list  of  reference  to  the  current  literature. 

The  aim  of  the  paper  is  to  provide  a  survey  of  a  part  of  this  new  field, 
with  emphasis  on  mathematics.  We  present  proofs  of  essential  theorems 
because  we  feel  these  proofs  are  essential  for  understanding  the  mathe¬ 
matical  aspects  of  meshless  methods,  which  has  approximation  theory  as 
a  major  ingredient.  As  always,  any  new  field  is  stimulated  by  and  related 
to  older  ideas.  This  will  be  visible  in  our  paper. 
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1  Introduction 

a.  A  Brief  Historical  Review  of  the  Numerical  Solution  of  Partial  Differential 
Equations. 

The  numerical  solution  of  partial  differential  equations  has  been  of  central 
importance  for  many  years.  Significant  progress  has  been  made  in  this  area, 
especially  in  the  last  30  years;  this  progress  is  directly  related  to  the  develop¬ 
ments  in  computer  technology.  Methods  such  as,  for  example,  the  finite  Element 
Method,  are  used  in  all  fields  of  application. 


2 


Although  significant  progress  has  been  made,  numerical  methods  for  the 
solution  of  differential  equations  are  still  often  based  on  heuristic  ideas,  and 
verified  by  numerical  experiments.  Mathematical  analysis  is  often  shallow,  and 
fails  to  fully  address  important  issues  that  arise  in  the  application  of  the  methods 
to  important  problems  in  engineering  and  science. 

There  are  three  classical  families  of  numerical  methods  for  solving  PDEs: 

1.  Finite  Difference  Methods 

2.  Finite  Volume  Methods 

3.  Finite  Element  Methods 

These  three  families  have  two  common,  basic  features: 

1.  They  employ  a  mesh; 

2.  They  use  local  approximation  by  polynomials. 

We  discuss  each  of  these  features  in  turn. 

Mesh  generation  is  often  very  expensive — especially  in  human  cost.  This  is 
for  several  reasons  for  this  cost: 

•  The  domain  of  the  underlying  problem  can  have  very  complex  geometry. 

•  The  domain  of  the  problem  may  change  with  time,  which  requires  remesh¬ 
ing  at  each  time  step,  as  for  example  in  the  problem  of  crack  propagation 
or  when  Lagrangian  coordinates  are  used. 

•  Adaptive  procedures  require  change  of  mesh  during  computation. 

Although  large  progress  has  been  made  in  the  theory  and  practice  of  mesh 
generation,  the  construction  of  the  mesh  is  still  a  very  delicate  component  of 
the  numerical  solution  of  differential  equations.  For  this  reason  there  is  an 
interest  in  the  development  of  methods  that  eliminate  or  reduce  the  need  for  a 
mesh. 

Although  polynomials  have  outstanding  approximation  properties,  there  are 
situations  in  which  they  are  not  effective.  We  mention  problems  whose  solutions 
are  not  smooth  in  the  sense  that  they  may  not  have  several  bounded  derivatives. 
For  such  problems,  there  are  sometimes  other  approximating  functions,  which  we 
will  refer  to  as  special,  that  are  effective.  The  classical  methods  are  not  flexible 
in  this  regard:  they  do  not  use  these  special  non  polynomial  approximating 
functions.  There  is  thus  an  interest  in  developing  and  analyzing  methods  that 
can  flexibly  use  these  special  approximation  functions. 

This  created  the  need  to  develop  methods  that  address  both  of  these  issues — 
the  elimination,  completely  or  partially,  of  the  need  for  mesh;  and  the  effective 
use  of  special  (non  polynomial)  special  approximating  functions.  The  inspiration 
for  such  methods  came  mainly  from  two  sources. 

The  first  of  these  sources  is  the  class  of  classical  particle  methods  that  arise 
in  physical  simulation  in  connection  with  the  Boltzmann  Equation  or  with  fluid 
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dynamics.  Particle  methods  attempt  to  describe  the  motion  of  the  atoms  or 
their  averages  (or  their  density)  in  Lagrangian  Coordinates  (see  [39],  [40],  [65], 
[66],  [68],  [69],  for  example). 

The  other  source  is  the  idea  of  interpolation  in  the  context  of  the  general 
Variational  Methods  (of  Galerkin  type).  With  these  methods  one  has  a  finite 
dimensional  space,  the  trial  space,  the  method  selects  an  approximation  from 
this  space,  and,  under  certain  general  conditions,  it  is  known  that  the  error  in  the 
approximation  is  no  larger  than  a  constant  times  the  error  in  best  approximation 
by  functions  in  the  trial  space.  Thus  the  quality  of  the  method  is  determined 
by  the  approximation  property  of  the  trial  space.  It  is  thus  natural  to  try  to 
find  a  trial  space  that  has  good  approximation  properties.  This  property  relates 
directly  to  interpolation  by  the  approximating  functions.  For  functions  in  one 
dimension  this  was  a  classical  issue  in  numerical  analysis,  and  starting  around 
1950  was  studied  in  higher  dimension  and  for  an  arbitrary  distribution  of  points. 
It  was  recognized  that  the  construction  of  trial  spaces  could  be  based  on  the 
idea  of  interpolation. 

b.  Meshless  Methods. 

Let  us  now  make  the  discussion  of  Variational  Methods  more  precise.  We 
consider  an  elliptic  PDF,  which  has  the  variational  or  weak  form, 

u  €  Hi;  B{u,v)  =  iF{v),  for  all  v  G  H2,  (1-1) 

where  Hi,H2  are  two  Hilbert  spaces,  B{u,v)  is  a  bounded  bilinear  form  on 
Hi  X  H2,  and  !F{v)  is  a  continuous  linear  functional  on  H2.  Under  certain 
general  condition  (the  inf-sup  or  BB  condition;  see  [5], [11]),  the  solution  u  is 
characterized  by  (1.1).  We  are  interested  in  approximating  u.  Toward  that 
end,  we  assume  we  have  two  finite  dimensional  space  Mi  C  Hi,  M2  C  H2  that 
satisfy  the  discrete  inf-sup  condition  (see  [11]).  The  approximate  solution  umi 
is  characterized  by 


umi  G  Ml,  B{umi,v)  =  iF{v),  for  all  v  G  M2.  (1.2) 


As  a  consequence  of  the  fact  that  Mi  and  M2  satisfy  the  discrete  inf-sup  con¬ 
dition,  we  know  that  the  approximation  umi  is  quasi-optimal,  i.e., 


||u  -  umiWhi  <  c  inf  ||u 

xeMi 


Xllffi- 


(1.3) 


We  note  that  there  are  delicate  problems  related  to  the  discrete  inf-sup 
condition  when  the  spaces  Mi  and  M2  are  different;  as,  e.g.,  in  the  case  of 
mixed  methods.  In  [3]  different  spaces  are  used  (without  mathematical  analysis 
of  the  discrete  inf-sup  condition) . 

We  thus  see  that  the  quality  of  the  approximation,  i.e.,  the  error  \\u—umi  Ijffi 
is  mainly  determined  by  the  approximation  property  of  the  trial  space  Mi,  i.e., 

by 


El 


inf 

xeMi 


||w 


Xllffi- 
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It  is  thus  natural  to  select  the  trial  space  Mi  so  that  Ei  is  small.  To  do  this 
effectively  one  should  use  whatever  information  is  available  on  the  solution  u. 
Note  that  with  a  general  variational  method,  as  we  have  formulated  it,  there  is 
no  mention  of  a  mesh.  Of  course,  one  may  use  a  mesh  to  construct  a  good  trial 
space;  that,  in  fact,  is  exactly  what  is  done  with  a  usual  finite  element  method. 
For  example,  the  trial  space  is  the  space  of  piecewise  linear  functions  over  a 
mesh. 

Meshless  Methods,  however,  either  avoid  the  use  of  a  mesh,  or  use  a  mesh 
only  minimally,  for  example,  only  for  the  numerical  integration.  The  Petrov- 
Galerkin  method  given  by  (1.2)  is  a  meshless  method  if  the  construction  of  Mi 
and  M2  either  does  not  require  a  mesh  or  requires  a  mesh  only  minimally.  Thus, 
in  designing  Meshless  Methods  within  the  framework  of  variational  methods,  we 
have  two  general  goals: 

1.  The  construction  of  trial  spaces  Mi  that  effectively  approximate  the  so¬ 
lution,  and  the  construction  of  test  spaces  M2  ensuring  inf-sup  (stability) 
condition. 

If  the  solution  has  special  features,  e.g.,  if  it  is  not  smooth,  we  should 
have  the  flexibility  to  use  special  approximating  functions. 

2.  The  minimizing  of  the  need  for  a  mesh. 

In  meshless  methods,  there  is  sometime  a  mesh  in  the  background, 
used  for  numerical  integration,  but  one  may  not  need  a  mesh  gener¬ 
ator. 

We  note  that  there  are  meshless  methods  that  are  not  of  the  type  given  by  (1.2), 
e.g.,  methods  based  on  collocation,  but  the  construction  of  approximating  space 
follows  the  guidelines  of  the  construction  of  Mi  as  mentioned  before. 

The  approximating  (trial)  spaces.  Mi,  can  be  the  spans  of  specific  approxi¬ 
mating  functions  (shape  functions),  with  either  global  or  local  supports.  Polyno¬ 
mials  and  radial  functions  are  examples  of  approximating  functions  with  global 
supports.  See  [63]  for  a  discussion  of  the  use  of  polynomials  and  [27],  [74]  for  a 
discussion  of  the  use  of  radial  functions.  Another  type  of  approximating  func¬ 
tions  is  related  to  interpolation  and  data  fitting  procedures.  For  a  survey  of 
various  approaches  we  refer  to  [3],  [31],  [37],  [38],  [41],  [53],  [54],  [61],  and  [77]. 
Typical  finite  element  approximating  functions  and  spline  functions  have  local 
supports.  In  [16]  shape  functions  that  are  effective  for  the  approximation  of 
solutions  of  elliptic  equations  with  rough  coefficient  were  identified  and  ana¬ 
lyzed;  the  idea  in  [16]  was  extended  and  developed  in  [18].  The  approximating 
functions  used  in  [16]  and  [17]  can  be  characterized  as  solutions  of  particular 
homogeneous  differential  equations.  In  one  dimension,  L-splines  -  a  general¬ 
ization  of  splines,  satisfy  a  differential  equation  and  are  used  as  approximating 
functions  ;  see,  e.g.  [87].  Principles  for  the  selection  of  shape  function  were 
addresses  in  [12]. 
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We  note  that  in  the  engineering  literature  many  names  are  used  for  methods 
that  differ  only  in  the  construction  of  shape  functions,  or  in  their  implemen¬ 
tation;  see,  e.g.,  [30]  and  others  ([84]).  For  a  survey  of  results  on  Meshless 
Methods  we  refer  to  [15],  [22],  [34],  [42],  [56],  [57],  and  [76]. 

One  of  the  major  problems  of  meshless  methods  is  the  imposition  of  bound¬ 
ary  conditions,  especially  the  Dirichlet  boundary  conditions.  It  is  well  known 
that  if  the  underlying  problem  is  a  Dirichlet  BVP,  the  essential  boundary  con¬ 
dition  is  addressed  with  a  method  such  as  the  penalty  method  or  the  method  of 
Lagrange  multipliers.  On  the  other  hand,  the  boundary  condition  of  a  Neumann 
Problem  is  natural  and  doesn’t  need  to  be  explicitly  imposed  in  the  variational 
formulation.  In  both  the  situations,  a  simple  uniform  mesh  on  a  rectangle  con¬ 
taining  the  domain  can  be  used;  the  mesh  need  not  conform  to  the  boundary 
and  a  mesh  generator  is  not  needed.  These  ideas  are  classical  and  have  been 
extensively  analyzed  (for  example,  see  [11]).  These  ideas  of  imposing  boundary 
conditions  can  be  used  in  the  context  of  meshless  methods  and  this  approach 
was  also  mentioned  in  [56].  The  boundary  of  the  domain  does  come  into  play 
in  the  construction  of  the  stiffness  matrix,  but  a  mesh  generator  is  not  needed. 
This  approach  was  generalized  and  used  together  with  the  ideas  in  [16],  [17], 
and  [83]  in  solving  problems  with  very  complex  geometries;  see,  e.g.,  [82]. 

We  mention  finally  a  meshless  method — The  Generalized  Finite  Element 
Method  (GFEM) — that  attempts  to  simultaneously  achieve  the  two  goals  of 
variationally  formulated  meshless  methods.  With  this  method  one  begins  with 
a  Partition  of  Unity.  Gonstruction  of  a  partition  of  unity  is  a  relatively  simple 
task.  It  can  be  done  by  various  means.  One  is  to  use  a  simple  mesh,  for  example, 
a  uniform  mesh,  and  use  the  associated  hat  functions  as  the  partition  of  unity. 
We  could  also  use  ideas  from  various  interpolation  procedures,  e.g.  the  Shepard 
Method.  It  is  essential  that  the  construction  can,  but  is  not  required  to,  utilize 
the  geometry  of  the  domain.  The  partition  of  unity  on  the  domain  is  obtained 
by  restriction.  The  partition  of  unity  functions  typically  have  compact  supports 
with  small  diameters. 

Then  we  multiply  the  partition  of  unity  functions  by  functions  that  are 
defined  separately  and  independently  on  the  supports  of  the  partition  functions. 
In  this  way  we  create  shape  functions  that  belong  to  and  can  be  used  in 

the  variational  method.  We  thereby  obtain  a  large  flexibility  in  the  construction 
of  shape  functions,  and  the  associated  trial  spaces.  This  fiexibility  can  be  used  to 
construct  approximations  that  utilize  the  available  information,  the  character  of 
a  singularity,  or  a  boundary  layer,  e.g.,  on  the  approximated  function  (solution). 
Hence  the  method  achieves  the  goals  mentioned  above. 

We  do  face  three  serious  difficulties  in  the  implementation  of  the  GFEM. 
First  there  is  the  problem  of  numerical  integrations  when  the  areas  over  which 
we  integrate  are  not  simple  triangles,  simplicies,  etc.,  as  with  the  usual  FEM.  We 
note,  however,  that  the  process  is  completely  parallelizable.  A  second  difficulty 
is  the  treatment  of  essential  boundary  conditions.  The  third  issue  concerns 
the  system  of  linear  equations.  It  may  be  singular,  and  thus  certain  classical 
methods,  such  as  multigrid,  may  not  be  applicable.  These  difficulties  can  be, 
and  have  been,  overcome  in  some  implementations,  so  it  is  clear  that  the  GFEM 
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shows  a  definite  advantage  over  the  classical  FEM  in  certain  situations.  We 
mention  problems  with  complex  geometry,  crack  propagation,  and  analysis  of 
mulit-site  local  damage. 

Of  course,  any  new  method  should  be  compared  with  previously  developed 
methods,  and  the  class  of  problems  for  which  the  new  method  is  superior  should 
be  identified.  Theoretical  and  practical  experience  (see  [15],  [56],  [83])  is  pro¬ 
gressing  in  this  direction.  Meshless  Methods  in  various  forms,  e.g.,  within  the 
framework  of  collocation  or  variational  methods,  are  now  the  subject  of  many 
papers  and  (engineering)  books,  which  mainly  focus  on  practical  aspects  without 
serious  theoretical  analysis. 

This  paper  focuses  on  ideas  and  theoretical  results.  Some  are  adjustments 
of  old  ideas  and  results.  Some  results  are  based  on  papers  that  are  submitted 
or  in  the  final  stage  of  preparation.  Although  we  focus  on  the  theory,  we  have 
attempted  to  address  theoretical  issues  that  illuminate  practical  issues.  We  will 
show  that  the  results  presented  here  are  natural  generalizations  of  the  classical 
FEM,  which  is  a  special  case  of  some  of  the  methods  presented  here.  This  paper 
addresses  only  problems  related  to  linear  PDEs. 

Various  relevant  and  typical  references  are  provided.  The  reference  list  is  not 
comprehensive,  although,  together  with  the  citations  in  the  references  provide, 
in  our  opinion,  a  very  reasonable  description  of  the  current  state  of  the  art  for 
meshless  methods. 

c.  The  Scope  of  this  Paper. 

The  short  Section  2  defines  the  model  problem,  a  linear  elliptic  boundary 
value  problem.  Section  3.1  presents  approximation  results  when  the  particles 
are  uniformly  distributed.  The  presented  results  were  obtained  using  the  Fourier 
Transform.  Section  3.2  presents  an  alternate  proof  of  the  approximation  results 
that  can  be  generalized  to  the  case  of  non-uniformly  distributed  particles.  Sec¬ 
tion  3.3  discusses  approximation  for  arbitrarily  distributed  particles.  Section 
4  discusses  the  construction  of  shape  functions,  and  presents  some  results  on 
interpolation  and  on  the  asymptotic  form  of  the  error.  Section  5  addresses  the 
question  of  superconvergence.  Section  6  discusses  the  Generalized  Finite  Ele¬ 
ment  Method.  Section  7  discusses  the  application  of  the  approximation  results 
developed  in  Sections  3,  and  discusses  the  treatment  of  Dirichlet  boundary  con¬ 
ditions.  Section  8  explains  some  implementational  aspects.  Section  9  reports 
some  numerical  examples  obtained  by  the  GFEM,  when  the  domain  is  very 
complex.  Finally,  Section  10  presents  additional  results  and  challenges. 

2  The  Model  Problem 

For  concreteness  and  simplicity  we  will  address  the  weak  solution  of  the  model 


problem 

— Au  +  u  =  f{x),  on  n  C  i?" 

(2.1) 

and 

du 

(2.2) 

—  =  0  on  ail 
on 
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or 


M  =  0  on  dQ,  (2.3) 

where  /  G  L2(fl)  is  given.  We  will  assume  that  is  a  Lipschitz  domain;  addi¬ 
tional  assumptions  on  dfl  will  given  as  needed. 

The  weak  solution  ug  G  H^{n)  {Hq{Q),  respectively)  satisfies 

B{uo,v)  =  B{v),  for  all  v  G  {v  G  i7g(n), respectively),  (2.4) 


where 

B(u,v)  =  /  {Vu -Vv  +  uv)  dx  !F{v)  =  /  fvdx.  (2.5) 

Jo.  Jo. 

The  energy  norm  of  ug  is  defined  by 

II^^oIIb  =  B{ug,UgY^"^  =  ||  Ug  || (Q)  .  (2.6) 

We  will  write  H  instead  of  H^{^)  or  i^g(^2)  if  no  misunderstanding  can  occur. 

Let  S'  C  iL  be  a  finite  dimensional  subspace,  called  the  approximation  space. 
Then  the  Galerkin  approximation,  ug  G  S,  to  ug  is  determined  by 

B{us,v)  =  B{v),  for  all  v  G  S,  (2-7) 

where  B  is  either  i?  or  a  perturbation  of  B.  If  B  =  B,  it  is  immediate  that 

ll^^o  —  W5||//i(n)  =  inf  ||mo  —  x||//i(o).  (2.8) 

Hence,  the  main  problem  is  the  approximation  of  ug  by  functions  in  S. 

Remark  2.1  The  Finite  Element  Method  (FEM)  is  the  Galerkin  Method  where 
S  is  the  span  of  functions  with  small  supports.  For  the  history  of  the  FEM,  see 
[9]  and  the  reference  therein. 

Remark  2.2  The  classical  Ritz  method  uses  spaces  of  polynomials  on  for 
the  approximation  spaces;  see,  e.g.,  [63]. 

As  mentioned  above,  the  Finite  Element  Method  uses  basis  functions  with 
small  supports,  e.g.,  “hill”  functions.  The  theory  of  approximation  with  general 
hill  functions  with  translation  invariant  supports  was  developed  in  1970  in  [4] 
using  the  Fourier  Transform.  The  results  in  [4]  were  applied  to  the  numercial 
solution  of  PDE  in  [5].  A  very  similar  theory,  also  based  on  the  Fourier  Trans¬ 
form,  was  later  developed  in  [80]  and  [81];  see  also  [55].  Later,  hill  functions 
were,  in  another  context,  called  Particle  Functions  (see  [39]).  In  the  1990s,  hill 
functions  began  to  be  used  in  the  framework  of  meshless  methods.  For  a  broad 
survey  of  meshless  methods  see  [56].  A  survey  of  the  approximation  properties 
of  radial  hill  functions  is  given  in  [27]. 

In  this  paper  we  will  survey  basic  meshless  approximation  results  and  their 
use  in  the  framework  of  Galerkin  Methods. 

8 


3  Approximation  by  Local  Functions  in  R^;  the 
/i-version  Analysis 

As  mentioned  in  Section  2,  we  are  interested  in  the  approximation  of  functions 
by  particle  shape  functions.  We  first  consider  uniformly  distributed  particles, 
and  then  general — non-uniformly  distributed — particles. 


3.1  Uniformly  Distributed  Particles  and  Associated  Par¬ 
ticle  Shape  Functions 

Let 

=  {j  =  (ji,  •  ■  •  ,jn)  integers} 

be  the  integer  lattice,  and,  for  j  =  (ji,  •  •  •  ,jn)  G  Z”  and  0  <  /i  <  1,  let 

=  {jih, . . .  ,jnh)  =  hj. 

The  points  are  called  uniformly  distributed  particles.  When  considering  such 
families  of  particles,  we  often  construct  associated  shape  functions  as  follows. 
Let  (j)  G  for  some  0  <  g,  be  a  function  with  compact  support;  let 

rj  =  supp  (j),  and  suppose 

ry  C  Bp  =  {x  G  M”  :  ||x|p  =  x\a - h  x^  <  p}. 

We  assume  that  0  G  ry  (iy  is  the  interior  of  ry).  Then  define 


=  (j)'^{xi,---  ,x„) 


/  x-jh\ 

V  h  J 


/xi  -  jih 

V  h 


for  j  G  Z"  and  0  <  h  <  1.  Clearly, 

T]^  =  supp  (j)'^  =  |x  :  ^  G  jyj  C  =  {x  :  ||x  -  x'ljW  <  ph}, 


and  x^  G  fj^.  Particle  and  particle  shape  functions  defined  in  this  way  are 
translation  invariant  in  the  sense  that 

x^j+i  =  x’]  +  xf  and  4>^j+i{x)  =  (/)^ (x  -  xj^) , 

and  will  sometimes  be  referred  to  as  translation  invariant.  They  are  a  special 
case  of  general  (non-uniformly  distributed)  particles,  which  will  be  addressed 
in  Section  3.3.  We  refer  to  ft.  as  the  size  of  the  particle  and  the  function  (j) 
is  called  the  basic  shape  function.  In  this  section  we  will  be  interested  in  the 
approximation  properties  of 
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which  is  the  linear  span  of  the  associated  shape  functions,  as  — >  0.  The 
parameter  k  in  is  related  to  a  property  of  (/)^(a;)’s,  which  will  be  discussed 
later.  We  will  refer  to  as  the  particle  space  in  M”.  The  w^' s  are  called 
weights.  Specifically,  given  u  G  we  are  interested  in  estimating 

inf  IIu-xIIhoCK"),  (3.3) 

for  0  <  s  <  min{g,  k  +  1}.  We  are  especially  interested  in  the  maximum  fj,  such 
that 

inf  ||u- x||^s(R«)  <  C'(fc,g)/i^||'u||jjf.+i(R„),  (3.4) 

for  0  <  s  <  minj^,  k  +  1},  where  the  constant  C  =  C{k,  q)  depends  on  k,  q,  but 
is  independent  of  ft-  (C  also  depends  on  cj)). 

Because  we  are  assuming  the  particles  are  uniformly  distributed,  and  hence 
the  particles  and  shape  functions  are  translation  invariant,  estimates  of  the  form 
(3.4)  can  be  obtained  via  the  Fourier  Transform.  This  was  done  in  [4]  and  [80], 
[81].  We  will  cite  one  of  the  results  in  [81]. 

Let 

<^(C)  =  0(?ir--  ,Cn)  =  /  </)(x)e"“'^ da; 

JM" 

denote  the  Fourier  Transform  of  We  note  that  G  (^“(IR”)  since  (f>{x) 

has  compact  support.  We  use  the  usual  multi-index  notation:  a  =  (oi, . . . ,  a„), 
with  at  >  0;  |q;|  =  ai  -I-  •  •  •  -I-  a„;  a;“  =  a;“^  •  •  •  a;“";  and 

- , 

^  aer  •  •  • 

Theorem  3.1  [81]  Suppose  <j)  G  iL*(IR”)  has  compact  support,  where  the  smooth¬ 
ness  index  q  >  0  is  an  integer.  Then  the  following  three  conditions  are  are 
equivalent: 

1. 

m  +  0  (3.5) 

and 

D°‘(j>(2'Kj)  =  0,  for  0  yf  j  G  Z"  and  |q;|  <  k,  (3.6) 

where  k  is  a  non-negative  integer. 

2.  For  |q;|  <  k, 

j°‘4>{x  —  j)  =  dx“  -I-  q°‘{x),  for  all  x  G  M”,  (3.7) 

jeZ" 

where  d  yf  0  and  q°‘{x)  is  a  polynomial  with  degree  <  |a|. 

The  equality  in  (3.7)  is  equality  in  L2(IR”);  i-e.,  equality  for  almost  all  x  G 
M".  The  function  of  the  right-hand  side  of  (3.7)  is,  of  course,  continuous. 
If  the  function  on  the  left-hand  side  is  continuous,  which  will  be  the  case 
if  q  >  n/2,  then  (3.7)  will  hold  for  all  x  G  M". 
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3.  For  each  u  €  there  are  weights  Wj  G  M,  for  j  G  Z"  and  0  <  h, 

such  that 

iim  - 

<  /or  0  <  s  <  min{g,fc  +  1},  (3.8) 

and 

(3-9) 

jeZ" 

Flere  C  and  K  may  depend  on  q,  k,  and  s,  but  are  independent  of  u  and  h.  The 
exponent  fc  +  1  —  s  is  the  best  possible  if  k  is  the  largest  integer  for  which  (3.7) 
holds. 

If  (3.7)  holds,  the  basic  shape  function  (j)  is  called  Quasi- Reproducing  of 
Order  k.  If  (3.7)  holds  with  d  =  1  and  q°‘{x)  =  0,  ^  is  called  Reproducing  of 
Order  k.  If  </>  is  quasi-reproducing  of  order  k  (respectively,  reproducing  of  order 
k),  then  the  corresponding  particle  shape  functions  (jf)  are  also  called  quasi- 
reproducing  of  order  k  (respectively,  reproducing  of  order  k).  The  parameter 
k  in  defined  in  (3.2),  is  the  quasi-reproducing  order  of  the  basic  shape 

function  </>. 


Remark  3.1  If  one  were  to  define  the  notion  of  quasi-reproducing  basic  shape 
function  </>  of  order  k  with  the  formula 

j°‘(t){x  —  j)  =  daX°‘  +  q°'{x),  for  all  x  G  M",  for  |a|  <  k,  (3.10) 

jGZ" 

where  da  yf  0,  it  might  appear  that  this  would  lead  to  a  larger  class  of  (j>’s.  This, 
however,  is  not  the  case;  it  is  easily  shown  that  if  (j)  satisfies  (3.10),  then  da  =  d, 
ior  \a\  <  k. 

In  one  dimension  we  can  prove  more. 


Theorem  3.2  [81]  Suppose  (p  G  has  compact  support  and  satisfies  Con¬ 

dition  1  in  Theorem  3.1,  i.e.,  satisfies  (3.5)  and  (3.6).  Then 


m 


m) 


|^sin(C/2)^"+^ 


(3.11) 


where  is  an  entire  function. 


Proof.  Because  p  has  compact  support,  ([{f,)  is  an  entire  function,  and 
because  of  (3.5)  and  (3.6),  (/(O)  yf  0  and  0(/)  has  zeros  of  at  least  order  k  at 
27r/,  0  yf  /  G  Z.  Let 


^fc(C) 


|^sin(//2)y+^ 


(3.12) 
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The  function  is  entire  with  only  zeros  of  order  /c  +  l  at  27rj,  for  0  yf  j  G  -E. 

Hence 


z{o  =  ko/^kio 


is  entire,  and 


as  desired.  □ 


z{0 


|^sin(C/2)^^+^ 


(3.13) 


Theorem  3.3  [4]  Suppose  (j)  G  has  compact  support  and  satisfies  Con¬ 

dition  1  in  Theorem  3.1,  i.e.,  satisfies  (3.5)  and  (3.6).  Then,  for  any  e  >  0, 


supp  (j)  (f 


(fc+1)  ,  ,  (fc+1) 
2  ’  2 


(3.14) 


Proof.  Suppose,  on  the  contrary,  that 


supp  (j)  C  [— (fc  +  l)/2  +  e,  {k  +  l)/2  —  e],  for  some  e  >  0.  (3.15) 


We  will  show  that  this  assumption  leads  to  a  contradiction. 

The  function  ^(^)  is  entire  and,  with  ^  +  z^2)  (3.15)  implies 


1^(01 


(3.16) 


This  estimate  follows  directly  from  the  definition  of  the  Fourier  Transform  and 
assumption  (3.15).  Using  elementary  properties  of  the  sine  function,  we  find 


that 


/  sin(g/2) 

V  e/2 


fc+1 


for  1^2!  large. 


(3.17) 


Using  (3.5),  (3.6),  (3.16),  and  (3.17),  we  have 


\zm 


<Co  +  Ck+i\^\'^+\  foralUeC, 


(3.18) 


where  Z{^)  is  as  in  (3.11).  Since  Z{^)  is  entire,  estimate  (3.18)  implies  (via  a 
generalization  of  Liouville’s  Theorem  for  entire  functions)  that  Z{()  is  a  poly¬ 
nomial  of  degree  <  fc  -I-  1.  Next,  we  use  (3.11)  and  (3.16)  to  get 


z{0 


|^sin(g/2)^^+^ 


|(/(C)|  < 


(3.19) 


Combining  this  estimate  with  the  lower  bound  in  (3.17)  we  have 
\Z{0\  <  for  1^2!  large. 


(3.20) 
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This  implies  Z{^)  =  0.  Thus,  (3.11)  implies  =  0,  which  contradicts  (3.5). 
Thus  (3.15)  is  false,  which  proves  (3.14).  □ 

The  case  of  uniformly  distribute  particles  is  very  special,  but  we  have  consid¬ 
ered  it,  and  cited  Theorem  3.1  from  [81]  because  the  result  provides  necessary 
and  sufficient  conditions  on  the  basic  shape  function  (j)  for  the  validity  of  the 
approximability  result  (3.8)  and  (3.9),  leading  to  the  optimal  value  for  /r  in 
(3.4). 

Comments  on  Theorems  3. 1-3. 3 

1.  ^(0)  yf  0  means  that  <j){x)  dx  yf  0. 

2.  For  the  validity  of  (3.8)  and  (3.9)  we  need  a  polynomial  reproducing  prop¬ 
erty,  namely  that  given  in  Condition  2.  Note  that  this  conditions  differs 
from  the  usual  polynomial  reproducing  property,  which  reads 

j°'4>{x  —  j)  =  a:“,  for  all  |q;|  <  p, 

ieZ" 


or,  equivalently, 

p{j)(j){x  —  j)  =  p{x),  for  all  polynomials  p{x)  of  degree  <  k. 


3.  Condition  2  implies  that 


<i){x  -  j)  =  d.  (3.21) 

Hence  the  functions  j  G  Z”,  are  a  partition  of  unity.  The  sets  fjj 

are  an  open  cover  of  M". 


4.  Taking  q  =  k  +  1  allows  application  of  Theorem  3.1  for  all  s  <  fc  -I-  1. 
Taking  q  >  k  +  1,  i.e.,  assuming  extra  smoothness  on  the  particle  shape 
functions  does  not  change  the  estimate. 


5.  Condition  (3.9)  is  a  stability  condition. 


6.  The  simplest  basic  shape  function  4>{x)  corresponds  to 


(3.22) 


Thus  they  (the  (/)’s)  are  spline  functions  over  a  rectangular  mesh,  which  are 
convolutions  of  the  characteristic  function  of  the  set  =  (—1/2, 1/2)". 


13 


7.  In  one  dimension  (n 
to 


and  they  are  B-splines.  Furthermore,  as  shown  in  Theorem  3.2,  if  4){x) 
satisfies  (3.5)  and  (3.6),  then  is  the  product  of  (Tfc(^)  and  a  suitable 
entire  function  Z{^).  Taking  account  of  the  growth  of  Z(^),  we  see  that 
it  is  the  Fourier  transform  of  a  function  with  compact  support.  Writing 
Z{^)  =  we  can  express  (3.11)  as 

Thus  any  (j){x)  that  satisfies  (3.5)  and  (3.6),  which  may  or  may  not  be 
piecewise  polynomial,  can  be  constructed  via  convolution  of  B-splines  with 
functions  of  compact  support.  If  n  >  1,  no  such  divisor  (j)/ak  exits  in 
general. 

8.  Theorem  3.3  has  as  especially  simple  interpretation  for  (/)’s  that  satisfies 

(3.5)  and  (3.6),  and  whose  support  is  an  interval.  Namely,  supp  ([)  D 
[  —  and  hence  grows  with  k.  As  is  well  known,  the  support 

of  the  B-spline  of  order  A:  is  [  —  and  hence  it  has  a  minimal 

support. 

9.  The  space  is  a  —regular  system  (this  notion  will  be  introduced 

in  Section  3.2),  with  k*  =  q  and  t  =  k  +  1.  -regular  systems  are 

analyzed  in  [11].  They  have  many  important  properties,  some  of  which 
will  be  used  in  the  following  sections. 

10.  The  approximability  of  the  classical  finite  element  shape  functions  (the 
hat  functions)  can  be  analyzed  with  Theorem  3.1  with  q  =  k  =  1. 

11.  The  weights  in  (3.8)  depend  on  u,  but  they  are  not  unique.  We  note  that 
the  functions  (f>j  may  be  linearly  dependent. 

3.2  Alternate  Proof  for  Uniformly  Distribnted  Particles 
and  Particle  Shape  Functions 

In  this  section  we  first  give  an  alternative  proof  that  Condition  2  in  Theorem  3.1 
implies  estimate  (3.8),  again  for  uniformly  distributed  particles  and  associated 
shape  functions.  This  alternative  proof  does  not  use  the  Fourier  Transform, 
and  it  can  be  naturally  generalized  to  the  non-uniformly  distributed  particles 
situation. 

We  review  our  notation  before  stating  the  theorem.  Recall  that 
x'j  =  jh,  for  j  G  Z"  and  0  <  h, 


=  1)  the  simplest  basic  shape  function  corresponds 


= 


sinU2 

U2 
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are  the  particles,  and  (j)  G  with  g  >  0,  is  the  basic  shape  function.  Also 

?7  =  supp  (f)  C  Bp,  and  0  G  fj.  Then  the  particle  shape  functions,  4>j{x),  are 
defined  by 

it  is  immediate  that 

r]f  =  supp  C 

and  x'j  G  fij. 

Theorem  3.4  Suppose  (p  G  with  smoothness  index  q  >0,  has  compact 

support  rj  C  Bp,  and  suppose  the  (f>j{x)  are  defined  in  (3.1).  Suppose  k  = 
0, 1, 2,  •  •  •  and  suppose,  for  |q;|  <  k, 

rfiix  -  j)  =  dx°^  +  q°‘{x)-,  (3.24) 

here  d  yf  0  and  q°‘{x)  is  a  polynomial  of  degree  <  \a\,  i.e.,  suppose  </>  is  quasi- 
reproducing  of  order  k.  Suppose  u  satisfies 

where  0  <rj  <  k,  (3.25) 

where  p  >  1  is  sufficiently  large  and  independent  of  h.  Then  there  exist  weights 
such  that 


lh\l2 
\H‘ 


<  C  for0<s<  min{g,rj  +  l}, 


jez„ 


ieZ"  j& 

where  C  is  independent  of  u  and  h.  If  u  G  +^(18"),  where  0  <  k'  <  k,  then 
\\u-  for0<s<  mm{q,k'  +  1}. 


(3.26) 


ieZ" 


(3.27) 


Proof.  The  proof  is  in  several  steps. 

1.  Suppose  (j)  satisfies  (3.24),  and  write  g“(a;)  =  X]|7|<|a|-i Then 


jGZ" 


Y  (a;) 

ieZ" 

jeZ" 


a 

)  +9“ 

1  \h^ 

V ft,/  J 

da;“  +  /rl“l 

E 

^70  (l) 

|7|<|q|-1 


da;“+  'Y  ^'^^djax'^,  for  |q;|  <  fc.  (3.28) 

|7|<|a|-l 
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Equations  (3.24)  and  (3.28)  are,  in  fact,  equivalent;  (3.28)  could  be  viewed  as 
a  scaled  version  of  (3.24).  For  any  p  G  ,  there  is  a  uniquely  determined 
w  =  Wp^h  GV^  satisfying 

P{^)  =  ^  Wp^h(xj)4>j(x),  for  all  x  G  M".  (3.29) 

We  first  prove  the  existence  of  Wp^h,  and  begin  by  considering  the  monic  poly¬ 
nomials:  Pa  =  x°‘.  Suppose  |a|  =  0.  Then  from  (3.28)  we  have 

1  H  W{^^h{x^j)4>']{x), 

ieZ"  jeZ" 


where  ri;{i}_;i(a;)  =  1/d.  Next  suppose  |a|  =  1.  Using  (3.28)  again  we  have 


»"  =  E  )(»?)“<(»)  -  % 


jGZ" 


where 


a;“  hd^a 

T  “ 


Proceeding  in  this  way,  by  induction,  we  get  W{a:<»},;i(a:)  for  |q:|  <  k,  where 
W{x‘^}^h{x)  is  of  the  form 


1«{a;“},?i(a;)  =  CaaX"  -f  ^  60/3/1'“'  (3.30) 

|/3|<|a|-l 

where  eco,  =  d~^  and  Cap  are  expressions  in  d.ya,  I7I  <  |q:|.  For  p{x)  = 
I]|a|<fcCax“,  we  let  Wp^h{x)  =  Z]|o|<fc  CoW{a;op/j(x).  It  is  immediate  that 


P{x)  =  Y  '^P,h(Xj)(f)j(x), 

jGZ" 

which  establishes  the  existence  of  Wp^h{x).  One  can  show  that 

Wp^hix)  =  Y  [  Y  606/0/3/1'“'“'^'] x^.  (3.31) 

|/3|<fc  |/3|  +  l<|a|<fc,  a—0 


To  prove  the  uniqueness,  suppose  Wp^h{x)  =  0.  We  will  show  that  p{x)  = 
X]|a|<fc  Cax“  =  0.  Since  Wp^h{x)  =  0,  it  is  clear  from  (3.31)  that  the  coefficient  of 
xP  is  zero  for  \f}\  <  k,  from  which  we  can  deduce  that  Cq,  =  0,  |a|  <  k,  and  thus 
p(x)  =  0.  It  will  be  convenient  to  write  Wp^h{x)  =  A’^p.  Then  ■'p^^-pk  jg 
a  1-1,  onto  mapping  satisfying 

p{x)  =  Y^  (-d^p)(x^)(/)^(x),  for  all  x  G  M",  for  any  p  G  .  (3.32) 

jeZ" 
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We  define  A  =  when  h  =  1.  We  note  that  A  satisfies  (3.32)  with  h  =  1.  We 
also  have 

[{A^)-^w]{x)  =  all  a;  G  K",  for  any  w  G  V'^ .  (3.33) 

It  is  also  clear  from  the  construction  of  that  A^  :  ^  for  i  <  k. 

2.  Define  the  cells  ojj  and  loj: 

u)j  =  {x  :  \\x  -  j\\oo  =  max  \xi  -  ji\  <  p} 

and 

uj'^  =  {x  :  \\x  -  x^Woo  =  max  \xi  -  x’^  l  <  ph}. 

The  families  {ujj}j^z^  and  are  open  covers  of  M”  provided  p  >  1/2. 

Let 

A'J  =  {I  G  Z"  :  n  io^  +  0}, 

and  define 


It  is  immediate  that  one  can  select  M  and  p  such  that 

card  <M  (3.34) 

and 

n’;  C  (3.35) 

The  constants  M  and  p  are  independent  of  j  and  h,  but  do  depend  on  (p-, 
specifically  on  p. 

For  any  I  G  Z",  since  u  G  iL’’'''"^(i?^^),  it  is  well  known  ([25], [26], [29])  that 
there  is  a  polynomial  of  degree  <  k  such  that 

\W-p‘’''\\h^(b‘.A  <  for  0  <  s  <  n  +  1  <  fc+  1,  (3.36) 

where  C  is  independent  of  u,  h,  and  I,  but  does  depend  on  k  (p*’^  can,  in  fact, 
be  chosen  such  that  its  degree  <  ri).  Define  the  weights 

wt  =  (^V’")(4).  (3.37) 

Let  j  be  fixed.  We  will  work  with  the  polynomial  p^'^,  which  satisfies  (3.36) 
with  I  =  j,  as  well  as  the  polynomial  p*’^.  Using  (3.37),  we  find 

11'*^  ~  ^ 

ieZ"  ^ 

-  11^  ~ 

leA^ 

+  l(-4V’'*)(x/)  -  (^V’")(^f)l  UnH^i^py  (3.38) 

leA^ 
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We  now  estimate  the  two  terms  on  the  right  side  of  (3.38). 

3.  From  (3.32)  and  the  definition  of  we  have 

=  51  (-^V)(a;f)<))f(a;),  for  X  e  uj^, 
ieZ"  i,=A^ 

for  any  p  G  .  Using  this  formula  and  (3.36)  with  I  =  j,  we  obtain  the  estimate 

||u-  ^(^V-")(xf)0fll^.(.^^)  = 

leA^ 

<  (3.39) 


for  the  first  term  of  (3.38). 

A  scaling  argument  shows  that 

Thus 

^  |AV’'*(xf)  -  -4V’"(xf)l 

iga^ 

<  l-4V’'*(4)  -  AV’^(xf)|.  (3.40) 

ieA>; 

It  remains  to  estimate  the  right  side  of  this  inequality. 

For  I  G  A^,  LUi  C  and  hence,  using  (3.35),  wf  C  Also  wf  C 
Thus,  using  (3.36)  with  s  =  0,  we  have 


< 

~  ^\\h°{b^-^)  +  b  ~  p’''^\\h»{b’^.^ 

< 

+  Ch'^‘'^^\\u\\Hr-i  +  l(^gl_^y  (3.41) 

It  is  easily  shown  that  there  is  a  constant  C  such  that 

lklL~(u.'*)  <  for  any  w  G  V'^;  (3.42) 

C  is  independent  of  w,  h,  and  1.  Applying  (3.42)  to  w  =  we 

have 


|(AV’")(xf)-(AV’")(xf)| 

<  C'/r-"/2|jAV’'‘  -  (3.43) 
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For  any  p  G  ,  we  write  p{x)  =  p  ^  where  p  G  V^.  Using  (3.32)  with 
h  =  1  (recall  that  A  =  for  h  =  1),  we  see  that 


and  therefore 


ieZ" 


iGZ"  \  / 

ieZ" 

iGZ" 

iGZ"  ^  ^ 


Comparing  the  above  expression  with  (3.32)  and  using  the  uniqueness  of  the 
representation  (3.29),  we  obtain 


{A'^p){x)  =  (AP)  .  (3.44) 

We  further  note  that  is  a  norm  on  p,  and  since  all  norms  are  equiv¬ 

alent  on  a  finite  dimensional  space,  we  have 

W^pWhom  ^  C'IIpIIh»(‘^o)-  (3.45) 

Therefore  from  (3.44)  and  (3.45),  and  using  the  transformation  y  =  {x  —  x^)/h, 
we  have 


WA^pW 


2 


\{A’^p{x)\’^  dx 


\{Ap{y)\^dy 


<  C”  /  \p{y)fdy 

JUJQ 

=  c[  \p{{x- x'l)/h)\'^dx 
JA 


=  C  [  \p{x)\‘^  dx 
Juj^ 

=  for  P  e  (3.46) 
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with  C  independent  of  p,  I,  and  h.  Combining  (3.43),  (3.46)  with  —  p^'^, 
and  (3.41)  yields 

|(^V")(xf)-(-4V'")(xf)| 

and  hence,  using  (3.34),  we  have  the  estimate 

^  |(^V’")(xf)  -  (-4V’")(xf)l 

leA'j 

+  ^  (3.47) 

leA^ 

for  the  right  side  of  (3.38).  Now  we  combine  (3.38),  (3.39),  (3.40),  and  (3.47) 
to  obtain 

||«-  ^  wUiWHH^r)  ^  ^  E  (3.48) 

/ez"  leA^ 

4.  Finally,  we  estimate  ||m  —  X^/eZ" llwiR")-  Using  (3.48),  which  is 
valid  for  all  j  G  Z”,  and  (3.34)  we  obtain 

Ik-  J2 ^  E  ik“  E 

/6Z"  jeZ"  ieZ"  ^ 

ieZ" 

<  C  J2  (3.49) 

iez- 

where  C  is  independent  of  u  and  h.  This  proves  (3.26). 

Suppose  u  €  +^(IR"),  where  0  <  k'  <  k.  Then  taking  rj  =  k'  in  (3.49), 

and  using  the  fact  that  the  overlap  in  {73^^}j6Z"  is  bounded  independently  of 
h,  we  get 

||«-  f  E  Ikll 

/eZ"  \ieZ" 

<  Ch^  *lkllBfc'  +  i(K"); 

where  C  is  independent  of  u  and  h,  which  is  (3.27).  □ 

Remark  3.2  Estimates  (3.26)  and  (3.27)  have  been  established  provided  p  is 
sufficiently  large;  specifically,  provided  (3.34)  holds.  As  pointed  out  in  connec¬ 
tion  with  (3.35),  p  depends  on  p.  Note  that  the  constants  C  in  (3.26)  and  (3.27) 
depend  on  p. 
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So  far  in  this  section,  we  have  considered  functions  u  defined  on  M",  and  have 
presented  a  result  on  the  approximation  of  u  by  particle  shape  functions.  We 
now  consider  functions  u  defined  on  a  bounded  domain  in  M"  with  Lipschitz 
continuous  boundary.  We  will  show  that  defined  in  (3.2),  when  restricted 
to  Q.,  provides  accurate  approximation  to  u. 

We  first  recall  the  well  known  extension  result  ([78]),  that  there  is  a  bounded 
extension  operator  E  :  L2{fl)  — >  L2(M"),  i.e.,  an  operator  E  satisfying  Eu\n  =  u 
for  all  u  €  L2(fi),  such  that  if  m  €  i7'"(n),  then  Eu  €  i7™(M")  and 

\\Eu\\H^{Mr.)  <  Cm\\u\\H^{Q),  for  all  M  G  77™ (ff ) ,  w  =  0,  1,  . . . .  (3.51) 

Here  Cm  is  independent  of  u  but  depends  on  m. 

We  define  a  subset  of  Zq  of  Z",  which  will  be  used  in  the  next  result,  by 

=  {j  G  Z"  :  n  ^  0},  (3.52) 

where  r/j  =  supp 

Theorem  3.5  Suppose  (p  G  77'?(IR"),  with  smoothness  index  q  >0,  has  compact 
support  Tj  C  Bp,  and  is  quasi-reproducing  of  order  k.  Suppose  u  G  ~'’^(f^)> 
where  0  <  k'  <  k.  Then  there  are  weights  Wj  such  that 

\\u-  0  <  s  <  min(9,A:'  +  l),  (3.53) 

where  C  is  independent  of  u  and  h. 

Proof.  Suppose  u  G  ^^(f^))  and  let  u  =  Eu,  where  E  is  the  extension 
operator  mentioned  above.  Applying  (3.27)  of  Theorem  3.4  to  u,  there  are 
weights  such  that 

||m-  Y  ^  Ch^'+'^-^\\u\\Hk'+^^r,y 

Therefore,  from  (3.51)  with  m  =  k'  we  have 

\\u-  Y  =  11“-^ 

zeZ"  /eZ" 

<  11^  ~  T 

/eZ" 

<  Ch^  *l|ii||_f/'='+i(R") 

<  Ch^  (3.54) 

From  the  definition  of  Zq  in  (3.52),  it  is  clear  that 


\\u-Y^Ui\\HHn)  =  \\u 
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and  therefore  from  (3.54),  we  get  the  desired  result.  □ 

By  examining  the  approximation  of  u,  obtained  in  Theorem  3.5,  namely 

we  see  that  the  sum  involves  only  those  /’s  for  which 

supp  n  n  0, 

i.e.,  only  those  particles  x'l  such  that  dist(xf,n)  <  ph.  So  the  approximation 
involves  particle  shape  functions  corresponding  to  the  particles  inside  as 
well  as  some  particles  lying  outside  Q,.  We  will  denote  the  span  of  these  shape 
functions  by 

^n,h  =  span{<(i^|j.2  :  supp  0}.  (3.55) 

Thus  the  functions  in  are  the  functions  in  restricted  to  fl. 

(t, k*)-regular  systems: 

We  now  introduce  the  notion  of  a  (t,  fc*)-regular  system  of  functions  (c/. 
[11]).  Let  n  C  M",  and  suppose  Sh{^)^  0  <  <  1,  is  a  one-parameter  family 
of  linear  spaces  of  functions  on  fl.  For  0  <  k*  <  t,  Sh{^)  will  be  called  a 
(t,  k*)-regular  system  and  will  be  denoted  by  (fl)  if 

1.  S'*’''*(fl)  C  for  0  <  /i  <  1.  (3.56) 

2.  For  every  u  G  with  0  <  I,  there  is  a  function  gh  G  such  that 

\\u- ghlln^Q)  <C!h>^\\u\\Hi(<:i),  for  0  <  s  <  min{Z,  fc*},  (3.57) 

where  p  =  minjt  —  s,l  —  s}.  The  constant  C  is  independent  of  u  and  h. 

We  now  introduce  two  additional  notions. 

(LA)  A  (t,  fc*)-regular  system  (fl)  is  said  to  satisfy  a  local  assumption  if 
for  u  G  iL^(n),  with  support  S,  the  function  gh  G  (fl)  in  (3.57)  can 
be  chosen  so  that  the  support  Sh  of  gh  has  the  property  that 

ShCS^’^  =  {xGn:  d{x,  S)  <  Xh}, 

where  d{x,  S)  is  the  distance  from  x  to  S,  and  A  is  independent  of  h. 

(lA)  We  say  that  (fl)  satisfies  an  inverse  assumption  (c/.  [11])  if  there  is 
an  0  <  e  <  /c*  such  that 

Il5llff'-*(n)  <  for  all  k*-e  <  r  <  k*  and  all  g  G  5'*’'"  (fl), 

where  C  is  independent  of  h  and  g  (it  may  depend  on  k*  and  e). 
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A  {t,  fc*)-regular  system  is  known  as  {t,  A:)-regular  system  in  classical  literature 
([11]).  We  have  decided  to  use  k*  instead  of  k  in  this  paper  for  notational  clarity. 

The  approximation  space  Vq’^,  defined  in  (3.55),  is  a  (t,  fc*)-regular  system; 
more  precisely  we  have, 

Theorem  3.6  Suppose  0  <  g  <  fc  +  1,  and  suppose  (j)  G  iJ*(IR")  has  compact 
support  and  is  quasi-reproducing  of  order  k.  Then  is  a  {k 1,  q) -regular 
system. 

Proof.  We  show  that  is  a  {t,  A:*)-regular  system  with  t  =  k  1  and 
k*  =  q.  Since  (f  G  iJ'^(IR”),  it  is  clear  that  C  iJ*(n)  and  thus  (3.56) 
is  immediate  with  k*  =  q.  Next  we  show  that  (3.57)  is  satisfied.  Suppose 
u  G  i7^(M”)  with  /  >  0.  If  /  =  0,  (3.57)  is  trivial.  So,  suppose  1  <  L  Applying 
Theorem  3.5,  specifically  (3.53)  with  k'  =  min(fc  +  1,1)  —  1,  we  get 

for  0  <  s  <  min{g,  min{l,  fc  +  1}}  =  min{/,g}  (since  g  <  fc  +  1),  where  p,  = 
minj/c  +  1  —  s,l  —  s}.  This  is  (3.57),  with  g^  =  t  =  k-\-  1,  and 

k*  =  q.  D 

We  now  show  that  satisfies  the  local  assumption  (LA). 

Theorem  3.7  Suppose  (j)  G  i7'^(M"),  where  0  <  g  <  fc  +  1,  has  compact  support 
and  is  quasi-reproducing  of  order  k.  Then  satisfies  the  local  assumption 
(LA). 

Proof.  Suppose  u  G  77* (fl)  such  that  supp  u  =  S  C  ft.  Consider  the 
approximation  of  u,  obtained  in  Theorem  3.5,  namely 

9h= 

A  careful  study  of  the  proofs  of  Theorems  3.5  and  3.4,  and  considering  the  zero 
extension  of  u  outside  LI,  reveals  that  for  j  G  Zq, 

Wj  =  0,  if  and  only  if  n  S'  =  0. 

Now  for  j  G  Zq  such  that  Wj  yf  0,  we  know  that  pj  =  supp  (fj  C 
Therefore,  Sh  =  supp  g/j  =  {x  G  IR”  :  d{x,S)  <  (p  +  p)h},  and  so  we  can  take 
A  =  (p  +  p)  in  the  definition  of  LA.  For  small  h,  we  have  Sh  C  LI.  Hence 
satisfies  the  local  assumption  LA.  □ 

Remark  3.3  The  particle  space  is  {k-\- 1,  g)-regular  and  satisfies  the  local 
assumption,  LA,  for  LI  =  K". 
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We  note  that  the  particle  spaces  and  V'^’*  will  also  satisfy  the  inverse  as¬ 
sumption  lA,  if  additional  conditions  are  imposed  on  the  shape  functions  }• 
We  will  formulate  these  conditions  in  Section  3.3  in  the  context  of  non-uniformly 
distributed  particles;  the  corresponding  conditions  on  the  shape  functions  as¬ 
sociated  with  uniformly  distributed  particles  can  then  be  obtained  as  a  special 
case. 

3.3  Approximation  by  Particle  Shape  Functions  Associ¬ 
ated  with  Arbitrary  (Non-Uniformly  Distributed)  Par¬ 
ticles  in  M”.  The  /i- version 

In  this  section  we  will  generalize  the  major  part  of  Theorem  3.4. 

Suppose  is  a  family  of  countable  subsets  of  M”;  the  family  is  in¬ 

dexed  by  the  parameter  v,  which  varies  over  the  index  set  N .  The  points  in  X'^ 
are  called  particles,  and  will  be  denoted  by  x,  to  distinguish  them  from  general 
points  in  M".  If  it  is  necessary  to  underline  that  x  G  X‘',  we  will  write  x  =  x''. 
To  each  gX  G  X''  we  associate 

•  /iji..  =  hJJ,  a  positive  number; 

•  a,  bounded  domain  in  K"; 

•  =  (j)%,  a  function  in  with  g  >  0  and  with  =  supp 

assumed  compact. 

The  numbers  will  be  referred  to  as  the  sizes  of  the  particles  x,  and 

the  functions  (j)^„  are  called  the  particle  shape  functions.  For  a  given  v  G  N ,  let 

=  {X\  {hi, 0^1,  ■ 

Xi'^  will  be  referred  to  as  particle- shape  function  system  —  and  {A4'^}i/gAr  as  a 
family  of  particle-shape  function  systems.  This  nomenclature  is  similar  to  that 
used  in  FEM  when  we  speak  of  a  triangulation  and  a  family  of  triangulation. 

Regarding  the  particle-shape  function  system,  we  make  several  assumptions: 
Al.  For  each  ly, 

U  <  = 

xeX" 

i.e.,  for  each  v,  is  an  open  cover  of  K”. 

A2.  For  X  G  X'’' ,  let 

Sl  =  {yGX'' 

There  is  a  constant  k  <  oo,  which  may  depend  on  {Xi^}i,^N,  but  neither 
on  n,  nor  on  a;  €  X'^ ,  such  that 

card  Sf  <  k,  for  all  x  G  X''  and  all  v  G  N. 
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A3.  For  all  x  G  X'' ,  for  all  G  A^,  a;  G  ryjl  and  fj^  C  uj^. 

A4.  For  X  G  X'',  let 

U 

y^Q'i 

where 

Ql  =  {y&X^ 

There  is  a  0  <  p  <  oo,  which  may  depend  on  {AT^j^gAr,  but  is  independent 
of  X  and  v,  such  that 

C  for  all  X  G  X'' ,  for  all  ly  €  N, 

where  B^f^^  is  the  ball  of  radius  centered  at  x,  namely, 

=  {a;  G  K”  :  l|a;  -  xll  <  ph'^J. 

A5.  For  each  x  G  X'',  there  is  a  one-to-one  mapping  :  7^^  — >  such  that 

^  G  and  any  p  G  'P'^,  (3.59) 

y6Qx 

and 

I!-4^p||l2((.j-)  <  C\\p\\L2(^^.p  for  all  p  G  for  all  y  G  Q^,  for  all  x  G 

A6.  For  any  0  <  s  <  q, 

\\ry\\HHu.i)  <  for  all  y  G  Q^. 

The  constant  C  may  depend  on  {AT^},yg7V,  but  is  independent  of  y,a;,  and 

V. 

A7.  There  is  a  constant  C  such  that 

lkllL“(a;p  <  C{h'Q-^/‘^\\w\\L2(^^u^p  for  any  w  G  P^ 
where  C  is  independent  of  y  and  v. 

Remark  3.4  From  the  definitions  of  Q'^  and  S'^,  and  assumption  A3,  it  is  clear 
that  Q'^  C  Hence  from  the  assumption  A2,  it  is  immediate  that 

card  Q'^  <  k.  (3.60) 

We  could,  of  course,  have  stated  (3.60)  as  an  assumption,  but  have  chosen  to 
state  card  5^  <  /t  as  an  assumption  because,  generally,  is  easier  to  work  with 
than  We  also  note  that  assumptions  A1-A4  imply  that,  for  any  x  G  M", 

card  {a:  G  X''  :  x  G  <  n,  for  all  v  £  N.  (3.61) 
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Remark  3.5  We  note  that  the  left-hand  side  of  (3.59)  in  A5  is  defined  for  all 
X  G  M",  but  the  equality  holds  only  for  x  G  wjl. 

Remark  3.6  Note  that  assumptions  A1-A7  can  be  thought  of  as  assumptions 
on  for  each  v  &  N]  they  are  assumptions  on  in  that  they  are 

assumptions  on  Ai''  for  each  v  and  that  the  constants  in  the  assumptions  do 
not  depend  on  ly. 

Remark  3.7  Assumptions  A5  effectively  defines  the  notion  of  quasi-reproducing 
shape  functions  of  order  k.  Note  that  the  condition  is  local:  the  operator 
depends  on  x,  the  sum  is  taken  only  over  y  G  and  the  equation  holds  only 
for  X  G  The  shape  functions  (j)'^  are  said  to  be  reproducing  of  order  k  if 

p{y)(j)y{x)  =  p(x),  for  X  G  M",  and  any  p  G  . 

yex- 

If  the  shape  functions  <j)'^  are  reproducing  of  order  k,  then  it  is  immediate  that 
they  satisfy  A5  with  AJi  equal  to  the  identity  mapping  for  each  x. 

Remark  3.8  Assumption  A5  implies 

=  M”,  for  each  v. 

xGX‘' 


Remark  3.9  Consider  uniformly  distributed  particles,  Xj,  and  associated  par¬ 
ticle  shape  functions,  (j)’-,  as  defined  in  Section  3.2.  Then  with  x‘'  =  x’-,  h'^  =  h, 
4'x  =  4’j}  E^nd  =  LOj,  as  defined  in  the  proof  of  Theorem  3.4,  the  asso¬ 
ciated  particle  shape  function  system  satisfy  assumptions  A1-A7.  Note  that 
A’i  =  A^h  =  A^  satisfies  A5. 

^  X. 

3 

Suppose  is  a  family  of  particle  shape  function  systems,  satisfying 

A1-A7.  Define 


=  span  :  x  G  for  each  v  &  N.  (3.62) 

The  next  theorem  gives  an  approximation  error  estimate  when  a  function  u, 
defined  on  M",  is  approximated  by  a  function  in  i/  G  N. 


Theorem  3.8  Suppose  the  family  of  particle- shape  function  systems  {AT^},ygAr 
satisfies  A1-A7,  and  hf  <  1  for  x  G  ,  v  G  N .  Suppose 


<  oo,  where  rf  <  k,  for  all  x  G  X'' ,  for  all  v  G  N, 


(3.63) 

where  p  is  introduced  in  A^.  Also,  suppose  that  operators  A^,  introduced  in  A5, 
satisfy 


||(A^-  for  allpG  y  G  Q^,  (3.64) 
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for  all  X  G  X'' ,  v  G  N,  where  C  is  independent  of  x  and  v.  Then  there  are 
weights  w'f  G  M,  for  y  G  X'^ ,  for  all  v  G  N ,  such  that 


II w  -  X!  ^y^yll^i‘=(R") 

yex- 


< 


C 


I/ex- 


1/2 


(3.65) 


for  0  <  s  <  vai{q,ry  +  1  :  y  G  X'',v  G  N}.  The  constant  C  depends  on  the 
constants  in  Assumptions  A1-A7  and  on  (3.64),  ^ut  neither  on  u,  nor  on  v. 


Note:  If  (ff’s  are  reproducing  of  order  k,  then  (3.64)  is  trivially  satisfied  {cf. 
Remark  3.7).  Since  in  practice,  mainly  shape  functions,  that  are  reproducing  of 
order  k,  are  used,  we  have  not  included  (3.64)  in  the  set  of  basic  assumptions 
(A1~A7). 

Proof  .  The  proof  of  this  result  is  parallel  to  the  proof  of  Theorem  3.4. 

1.  The  sets  ujf.  play  the  role  of  the  sets  coj  in  the  proof  of  Theorem  3.4.  The 

sets  Qf,  flf,  and  the  mapping  Af  play  the  roles  of  the  sets  A^f 

and  the  mapping  Aj,  respectively,  in  the  proof  on  Theorem  3.4.  Assumptions 
A1“A7  state  the  properties  of  these  sets  and  the  mappings  we  will  need  in  this 
proof. 

2.  For  any  y  G  X‘',  since  u  G  it  is  well-known  that  there  is  a 

polynomial  pH’’'  =  p^'  of  degree  <  k,  such  that 


h-Pk  I 


<  C{hA  y- 


+  1-S| 


(3.66) 


for  0  <  s  <  -I-  1  <  A:  -I-  1,  where  C  is  independent  of  u,  v  and  y,  but  does 

depend  on  k  (  p^'''  can,  in  fact,  be  selected  so  that  its  degree  ^  r^).  Define  the 
weights 

w'^=  {Alpy-niv),  (3.67) 

where  Ay  is  the  operator  introduced  in  assumption  A5. 

Let  X  G  X’'  be  fixed.  We  will  work  with  the  polynomial  p-’",  which  satisfies 
(3.66)  with  y  =  X,  as  well  as  the  polynomial  p^’‘^ .  Using  (3.67)  we  find 


I|m-  E 

yex- 
<  ||u- 

<||u- 

+  Y.  \{Alp^n{y)  -  (-4^1^’'^) (y) I  (3.68) 

y&Qy 
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We  now  estimate  the  two  terms  on  the  right-hand  side  of  (3.68). 

3.  From  the  assumption  A5,  we  know  that 

{■Kp){y)4'l{x)  =  p{x),  for  X  G  and  any  p  G 

v&Q'i 

Using  this  formula  and  (3.66)  with  y  =  x,  we  obtain  the  estimate 
\\u- <  \\u-p^’’'\\hhu.a 

veQl 

<  (3.69) 

'  pK' 

for  the  first  term. 

Using  assumption  A6,  we  have 

<  C(/.p-^+”/^  for  all  y  G  Q^. 

Thus 

^  \{Aip^ni.y)  -  {Aipyni.y)\ 

!/6Qx 

<  E  {h''^-^+-/^\{AljAn{y)  -  {Alpyn{y)\-  (3.70) 
2/6 


It  remains  to  estimate  the  right-hand  side  of  this  inequality. 

For  y  G  Q^,  we  have  ojy  G  and  hence,  using  the  assumption  A4,  uty  G 

Also  LOy  G  Thus,  using  (3.66)  with  s  =  0,  we  have 


<  ,  +  V  (3.71) 


'H  ^  y 

'■  ph'i'  - 


H  y-  (B“ 


Now,  using  the  assumption  A7,  we  have 


{A:iAn{y)-{A;pyn{y)\ 

<  I  [(A^  -  App^’Ky)  I  +  I  -  py-nm  I 

<  ||(A^-  +  \\A';{p^’''-pyn\\L^(u.i^) 

<C'(/rp-"/2{||(A^-App^’b!//oK) 

+\\A;{^^'^-pyn\\H0i.^J.  (3.72) 
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Also,  using  the  assumption  A5  and  (3.71),  we  obtain 


.  +  {KY^-^\\u\\  .^^^.  }.  (3.73) 

''  ph"!  -  -n  ^  ph}^! 


Moreover,  from  (3.64),  we  have 

and  from  (3.66),  with  y  =  x,  and  recalling  that  <  1,  we  get 


(3.74) 


lb^’^ll//0K)  <  Ib^’"  -  +  Ib||^0(^.) 


<  c{K) 


ip\r^  +  l  I 


■W 1 1 //’-!_+ 1(32- 


<  c'lbl!3''i+i(32 

1  Pi'S 


(3.75) 


Combining  (3.72)-(3.75),  we  have 

\{Aip^n{y)-{A;pyn{y)\ 


1  pK 


1  pf'^ 


,}• 

(3.76) 


Then  we  combine  (3.70),  (3.76),  (3.60),  and  assumption  A2  (c/.  (3.60))  to  get 

^  \{Aip^ni.y)  -  (-4;p^’b(y)|lb^lk=(.^) 

v^Q'i 


!/6<3£ 


<  ^  E 


i’„+i  I 


\A\\ H’'A\By,y 


(3.77) 


which  is  an  estimate  for  the  second  term  in  (3.68).  Thus,  from  (3.68),  (3.69), 
and  (3.77),  we  have 


ib-E 

y6Qx 


y6Qx 


3  (32 


(3.78) 


4.  It  remains  to  estimate  jb  —  J^yex-'  ll-tf“(R")-  Using  (3.78),  which  is 
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valid  for  all  a;  G  X'',  and  assumptions  Al,  A2,  A4,  we  obtain 
11'*^  ~  ^  llff=(R")  < 

yGX- 

< 

< 

which  is  (3.65).  □ 

It  will  be  useful  to  state  estimate  (3.65)  in  Theorem  3.8  in  certain  alternate 
ways.  Given  a  family  of  particle-shape  function  systems  satisfying 

A1-A7,  define 

h‘'  =  sup  for  each  v.  (3.80) 

xGX-'  ^ 

With  this  definition,  from  (3.65)  we  have 

ll«-  E  E  (3-81) 

yex-  yeX’' 

Now,  if  =  k' ,  where  0  <  k'  <  k,  for  all  y  €  X'',  for  all  then  (3.81)  leads 
to  the  following  result. 

Theorem  3.9  Suppose  the  family  of  particle- shape  function  systems  {AT^j^gAr 
satisfies  A1-A7,  (3.64),  in  addition,  suppose  h'^  <  1,  for  all  v.  Suppose 
||M||//fc'+i(Rn)  <  oo,  where  Q  <k'  <k.  Then  there  are  weights  rcj)  G  M  such  that 

ll«  -  E  +,(«„),  (3.82) 

yex- 

for  0  <  s  <  min(q,  k'  -\- 1),  where  C  is  independent  of  u  and  v. 

We  note  that  if  h^^  <  h''^,  vi,V2  G  N,  then  we  would  view  as  a 
“refinement”  of 

There  is  yet  another  way  to  state  the  estimate  (3.82).  Let  0  <  ft-  <  1, 
and  suppose  {AJ'^j^gAr,  a  family  of  particle-shape  function  systems  satisfying 
A1~A7,  (3.64),  and  in  addition, 

h''  =  sup  h'(.u  <  ft,  for  each  v.  (3.83) 

We  can  now  think  oi  v  =  v{h)  as  determined  by  ft,  although,  of  course,  many 
particle-shape  function  systems  satisfy  (3.83).  We  can,  in  fact,  think  of  having 
a  one-to-one  correspondence  between  u  and  ft.  Thus  we  can  regard  ft  as  the 
parameter  and  write  a  family  of  particle-shape  function  systems  as 

{AJ^}o</i<i  =  {A^,{ft^,  u’f,  'fr^x^xf^}o<h<l 


E  ii«-  S 

xex-'yeQ^  ’’  ^ 

C  E  )’  (3-79) 
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instead  of  With  this  understanding  that  1/  =  ^{h),  the  estimate 

(3.82)  can  be  written  as 

||u-  ^  (3.84) 

yex'* 

We  are  naturally  interested  in  having  h  \  0,  and  hence  in  considering  ^(h)’s 
for  which  \  0.  More  specifically,  we  will  often  consider  a  sequence  hi  > 

h2  >  ■  ■  ■  \,  and  the  corresponding  sequence  of  particle  systems  At ■  ■  ■ , 
where  vi  =  i'{hi). 

We  remark  that  the  estimate  (3.65)  is  stronger  than  (3.82)  and  (3.84),  in 
that  (3.65)  uses  h'^„  instead  of  the  larger  h'' ,  and  (3.65)  allows  a  more  gen¬ 
eral  regularity  assumption  on  the  function  u.  The  viewpoint,  outlined  in  this 
paragraph,  is  similar  to  the  usual  view  of  meshes  in  FEM. 

For  a  given  family  of  particle-shape  function  systems  ,  we  defined 

the  space  V*’'^  in  (3.62).  With  h,  0  <  /i  <  1,  as  the  parameter,  i.e.,  for  a  given 
family  0  <  h  <  1,  we  will  use  the  space 

^  =  span  :  X  G  (3.85) 

So  far,  we  have  discussed  the  approximation  of  a  function  u  defined  on 
M",  by  particle  shape  functions.  We  now  consider  u  defined  on  fl,  where  is 
bounded  domain,  with  Lipschitz  continuous  boundary,  in  K”.  We  now  show 
that  functions  in  defined  by 

:  (ji'l  G  where  ^  n  yf  0},  (3.86) 

provide  accurate  approximation  of  functions  u,  defined  on  fl. 


Theorem  3.10  Suppose  A4^,  0  <  h  <  1,  is  a  family  of  particle  shape  function 
systems  satisfying  A1-A7  and  (3.64-).  Let  Q.  C  IR”  he  a  hounded  domain  with 
Lipschitz  continuous  boundary,  and  suppose  u  G  +^(n),  where  0  <  k'  <  k. 
Then  there  are  weights  G  K  such  that 

ll«-  E  (3-87) 

yGX>' 


for  0  <  s  <  min(q,  k'  +  1),  where  the  constant  C  is  independent  of  u  and  h. 


The  proof  of  this  theorem  is  based  on  using  (3.84)  on  the  extension  u  =  Eu, 
and  is  similar  to  the  proof  of  Theorem  3.5.  We  omit  the  proof  of  this  theorem. 
We  note  that  the  approximation  obtained  in  Theorem  3.10,  is 

such  that 


Ehih 

WyCfy 

yex'* 


G  V 


k,q 

Q,h' 


In  Section  3.2,  we  reviewed  the  notion  of  (t,  fc*)-regular  system  Sh{Lt).  In 
the  next  theorem,  we  show  that  'Vq\  is  a  (fc  -I-  1,  g)-system. 
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Theorem  3.11  Suppose  0  <  h  <  1,  is  a  family  of  particle  shape  function 
systems  satisfying  A1-A7  and  (3.64)-  Let  LI  C  M”  be  a  hounded  domain  with 
Lipschitz  continuous  boundary.  Then  is  a  {k+  I,  q) -regular  system,  where 
k  is  the  order  of  quasi-reproducing  shape  functions  in  A4^. 

The  proof  of  this  theorem  is  similar  to  the  proof  of  Theorem  3.6,  and  will 
be  omitted. 

Remark  3.10  The  space  'Vq\  satisfies  the  local  assumption  (LA). 

Quasi-Uniform  Particle- Shape  function  System: 

We  will  call  a  family  of  particle-shape  function  systems  {Af^}o</i<i  quasi¬ 
uniform  if  there  is  a  /3,  1  <  /3  <  oo,  such  that 

h^ 

(3~^  <  ^  <  (3,  for  all  x,y  &  X^,  for  all  0  <  /i  <  1.  (3.88) 

hy 

We  note  that  (3.88)  is  equivalent  to 

for  all  y  G  X^,  for  all  0  <  /i  <  1,  (3.89) 

hy 

where  h  is  defined  by  (3.83). 

Remark  3.11  We  can  also  define  uniform  particle-shape  function  system  by 
imposing  the  condition 

h^  =  hy,  for  all  x,y  &  X^,  for  all  0  <  ft-  <  1. 

We  note  that  the  system  with  uniformly  distributed  particles  and  the  associated 
shape  functions  as  defined  in  Section  3.1  is  uniform.  But  uniform  particle  shape 
function  systems  may  have  particles  that  are  not  uniformly  distributed. 

Consider  a  family  of  particle-shape  function  systems  {Al^}o<h<i  satisfying 
the  assumptions  A1-A7.  Let  Ll  C  M"  be  a  bounded  domain,  and  define  A^  = 
{x  G  X^  :  fj^  n  n}.  Suppose  satisfies  the  following  additional  assumptions: 

•  A4^  is  quasi-uniform,  i.e.,  (3.89)  holds. 

•  For  all  X  G  A^,  there  is  a  /3  >  0  such  that,  for  0  <  s  <  g, 

<  f3h^~",  for  all  y  G  Q|,  (3.90) 

where  Q’f  =  {y  e  C  yf  0}  (c/.  A4). 
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•  For  all  Wy  G  M,  for  y  G  Q^,  and  x  G  A^,  there  is  C  >  0,  independent  of  x, 
such  that 

H  ^  ^*11  H  for  0  <  s  <  g.  (3.91) 

yeQj  y^^Q'i 

Then  the  particle  space  VqI  satisfies  the  inverse  assumption  lA,  introduced  in 
Section  3.2.  To  see  this,  consider  x  G  A^.  Then,  using  (3.90)  and  (3.91),  we 
have 

II  ^  ^  llr/9(whnn) 

yeQS  yeQS 

yeQ^ 

yeQ^ 

<  Ch^-^ll  ^  (3.92) 

yeQi 

where  C  depends  on  k  {cf.  A2).  Thus 

II  E  <  E  II  E 

y6^n  2/6  Q| 

<  E  II  E  ^yJyWln^inn) 

y&Q’i 

<  Ch^-n  Y. 

Since  any  element  g  of  Vq'(,  is  of  the  form  '^y^j^h  Wy(f>y\^,  we  have  shown  that 

'^n  \  satisfies  the  inverse  assumption  lA.  We  summarize  the  above  discussion  in 
the  following  theorem: 

Theorem  3.12  Suppose  ,  0  <  ft.  <  1,  zs  a  family  of  quasi-uniform  particle- 
shape  function  systems  satisfying  A1-A7,  (3.90),  and  (3.91).  Let  LI  C  M"  he 
a  bounded  domain  with  Lipschitz  continuous  boundary.  Then  Vq’^  satisfies  the 
inverse  assumption  I  A. 

Remark  3.12  We  can  show  that  the  particle  space  also  satisfies  the  inverse 
assumption  lA,  if  Ai^ ,  satisfying  A1-A7,  also  satisfies  (3.90)  and  (3.91)  with 
LI  =  M". 
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4  Construction  and  Selection  of  Particle  Shape 
Functions 

In  Section  3,  we  presented  an  abstract  description  of  particle-shape  function 
systems  with  respect  to  uniform  as  well  as  non-uniform  distribution  of  parti¬ 
cles.  We  showed  that  if  these  particle-shape  function  systems  satisfy  certain 
properties  (assumptions  A1-A7  and  (3.64)),  they  will  have  good  approximation 
properties.  In  this  section  we  will  present  an  example  of  a  particular  particle- 
shape  function  system,  where  the  shape  functions  are  reproducing  of  order  k, 
and  show  that  under  certain  conditions  they  satisfy  assumptions  A1--A7,  and 
hence  have  good  approximation  properties.  We  note  that  (3.64)  is  trivially  satis¬ 
fied.  This  example  will  also  show  that  a  wide  variety  of  particle  shape  functions 
can  be  constructed.  Therefore  it  is  important  to  address  the  issue  of  selecting  an 
appropriate  class  of  shape  functions  that  would  yield  efficient  approximation  of 
the  solution  of  a  particular  problem,  or  a  class  of  problems.  We  also  present  an 
interpolation  result  that  will  indicate  a  procedure  for  choosing  a  class  of  shape 
functions,  among  a  given  collection  of  such  classes.  Such  shape  functions  will 
yield  the  smallest  value  of  the  usual  Sobolev  norm  interpolation  error,  when  the 
interpolated  function  is  included  in  a  higher  order  Sobolev  space. 

4.1  An  Example  of  a  Class  of  Particle  Shape  Functions 

Several  particle  shape  functions  have  been  developed  over  the  last  decade.  SPH 
shape  functions  [39]  were  introduced  in  the  context  of  fluid  dynamics,  whereas 
Shepard  functions  [77]  and  MLS  shape  functions  [53]  were  introduced  in  the 
context  of  data  fitting  with  respect  to  irregularly  distributed  particles  in  higher 
dimensions.  In  the  90’s,  RKP  shape  functions  were  introduced  [58]  in  the  context 
of  approximation  of  solutions  of  partial  differential  equations.  In  this  paper,  we 
describe  the  construction  of  RKP  shape  functions  for  non-uniform  as  well  as 
uniform  distribution  of  particles,  and  relate  them  to  the  abstract  setting  given 
in  Section  3.  Specifically,  we  will  show  that  the  resulting  particle-shape  function 
system  satisfies  assumptions  A1-A7. 

N on-uniformly  distributed  particles: 

For  v  G  N,  N  an  index  set,  let  =  {xf }igz  where  G  M".  With 
each  xj'  €  X‘' ,  we  associate  a  positive  number  h^.  We  consider  a  fixed  i:  and 
often  suppress  the  superscript  v,  e.g.,  we  write  Xi  and  hi  instead  of  and  h'', 
respectively.  We  will  comment  about  v  when  appropriate. 

Let  w{x)  >  0  be  a  continuous  function  with  compact  support,  specifically, 

■q  =  supp  w{x)  =  BjiiQ),  i?  >  0.  (4.1) 

The  function  w{x)  is  called  a  weight  function  (or  window  function). 

The  commonly  used  weight  functions  in  1-d  are  as  follows: 
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(a)  Gaussian: 


\x\  <  R 
\x\  >  R, 


w{x)  = 


where  5  >  0. 

(b)  Cubic  spline: 


^5{x/R)^  _  gi5 


l-e-5  ’ 

0, 


(4.2) 


w{x)  = 


I  -  4,{x/RY  +  4(a;/i?)^, 
I  —  A{x/ R)  +  4(x/i?)^  — 

0, 


|x|  <  R/2 
i?/2  <\x\<R 
|x|  >  R. 


(c)  Conical: 


w{x) 


[l-{x/Rn,  \x\<R 

0,  |a:|  >  R, 


(4.3) 


(4.4) 


where  I  =  1,2....  We  note  that  one  may  consider  non-symmetric  versions  of 
some  of  these  weight  functions,  as  was  done  in  [2] . 

In  M",  w{x)  can  be  constructed  from  1-d  weight  function  w{x)  (symmetric) 
as  w{x)  =  wdlxll),  where  ||a;||  is  the  Euclidean  length  of  x.  w{x)  can  also  be 
constructed  as,  w{x)  =  nj=i'^G^)>  where  x  =  (ix,  2a:, . . . ,  „x)  G  M".  Conse¬ 
quently,  r]  will  be  an  n— cube.  However,  we  will  assume  rj  given  by  (4.1)  in  this 
section. 

For  each  j,  we  define 


Clearly, 


Let 


and  assume  that 


r]j  =  supp  Wj{x)  =  BRhjixj). 


Qi  —  ^Xj  .  T]i  yf  0}, 


(4.5) 

(4.6) 

(4.7) 


card  Qi  <  k,  for  all  f  G  Z, 


(4.8) 

(4.9) 


where  k  is  independent  of  i  and  v. 

For  a  given  integer  k,  k  >  0,  the  RKP  shape  function  4'j{x),  associated  with 
the  particle  Xj,  is  defined  by 

(j)j{x)  =  Wj{x)  ^  (x  -  Xj)°‘ba{x),  (4.10) 

\a\<k 

where  ba{x)  are  chosen  so  that 

''^p{xj)(j>j{x)  =  p{x),  for  X  G  ffi",  for  p  G  7^^(M"),  (4-11) 
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so  that  are  reproducing  of  order  k.  This  gives  rise  to  a  linear  system 

in  ba{x);  namely 


^  ma+fs{x)ba{x)  =  for  \f3\  <  k,  (4.12) 

\a\<k 

where  5\f3\,\a\  is  the  Kronecker  delta,  and 

'nia{x)  =  Wj{x){x  —  Xj)^. 
jez 

It  is  clear  from  (4.6)  and  (4.10)  that 

supp  (()j(a;)  =  swpp  Wj{x)  =  rij.  (4-13) 

We  now  briefly,  describe  the  derivation  of  (4.12).  For  a  fixed  y  G  M",  consider 
Pfi(x)  =  {y  -  x)^ ,  0<\p\<k. 

Using  p{x)  =  ppix)  in  (4.11)  we  get 

J2iy  -  Xjf(l)j{x)  =  {y-  xY, 

jez 

and  letting  y  =  a;  in  the  above  equality,  we  have 

^(a;-Xj)^^j(x)  =  (5|;3|,o,  0  <  |/3|  <  fc.  (4.14) 

Thus  (4.11)  implies  (4.14).  In  fact,  one  can  also  show  that  (4.14)  implies  (4.11). 
Now  using  (4.10)  in  (4.14),  we  get 


'^1/31.0  = 

lez 

=  “  Xjfwjix)  XI  (2^  “  XjYbYx) 

I  a  I  <  fc 

=  X  bYx)J2^jY)ix  - 

I  a  I  <  fc 

=  X  rna+p{x)baix),  (4-15) 

|a|<fc 

which  is  (4.12). 

We  now  consider  the  unique  solvability  of  (4.11).  For  k  =  0,  the  linear  system 
(4.12)  is  mo{x)bo{x)  =  [J2iez'^dY]bo{x)  =  1.  Assuming  Y.i&'^ii.x)  yf  0,  a;  G 
M",  we  have  bo(x)  =  l/mo(x).  Therefore  from  (4.10),  we  have 


^j(x) 


Wj(x) 

E^ez^i(Y’ 


j  G  Z. 
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This  expression  for  {(j)j(x)}  gives  another  verification  that  they  form  a  parti¬ 
tion  of  unity.  These  shape  function  are  called  Shepard  functions;  they  were 
introduced  in  [77]. 

The  unique  solvability  of  (4.12),  for  A:  >  1,  depends  on  the  weight  functions 
Wj’s  and  on  the  distribution  of  the  particles  {xi\  in  M".  The  required  distribu¬ 
tion  of  particles  in  turn  is  related  to  the  interpolation  problem  in  M".  It  was 
shown  in  [48]  that  a  necessary  condition  for  unique  solvability  of  (4.12)  is  that, 
for  X  e  K”, 

card  A{x)  >  dim  7^*,  (4.16) 

where 

A{x)  =  {xr.  x  &  rji}.  (4.17) 

For  fc  =  1,  in  was  shown  in  [48]  that  the  linear  system  (4.12)  is  non-singular  if 
the  following  conditions  are  satisfied: 

(a)  There  are  constants  Ci,  C2  >0  independent  of  v,  and  h  >  0,  such  that 

Cl  <  <  C2,  for  all  i  G  Z;  (4.18) 

h 

(b)  There  are  constants  Ci,C2  >  0,  independent  of  v,  such  that  for  any 
X  G  M",  there  are  (n  -I-  1)  particles  a;,,  G  A{x),l  =  0, . . . ,  n,  such  that 

min  w  (  ^  >  Cl  >  0  (4.19) 

o<;<n  \  h  J 

and 

Volume  K{xig,Xi^, . . .  ,Xi^)  >  C2h^,  (4.20) 

where  K{xi„,Xi^, . . .  ,Xi^)  is  the  simplex  with  vertices  xi, , I  =  0, 1, . . . , n. 

We  will  now  cast  RKP  shape  functions,  discussed  above,  in  the  framework 
of  a  particle-shape  function  system,  introduced  in  Section  3.3.  We  started  with 
a  collection  of  particles  X''  =  {xj}j^z,  where  Xj  G  M”,  and  positive  numbers 
hj.  Corresponding  to  each  particle  Xj  G  X‘',  we  associated,  in  (4.10),  the  RKP 
shape  function,  =  (j>j  with  compact  support  lyj  =  rjj  =  Buhjixj),  where 
the  parameter  R  was  related  to  the  weight  function  w{x).  It  was  shown  in 
[48]  that  if  w{x)  G  C*(M"),  then  (j)j  G  C'^(M"),  and  thus  (j)j  G  i7*(M”);  here 
we  assume  <7=1.  We  recall  that  the  conditions  (4.8),  (4.8),  (4.16),  (4.18)- 
(4.20)  were  required  for  the  construction  of  shape  functions,  (j>j,  j  G  Z.  We  let 
ujj  =  LOj  =  f/j;  certainly  wJ’s  are  bounded  domains.  We  now  show  that  the 
family  of  particle-shape  function  systems  {A4‘'}^(z]y,  where 

with  these  choices  for  (f)^  and  wf,  satisfies  assumptions  A1~A7  in  Section  3.3. 
We  will  continue  to  use  the  notation,  introduced  earlier  in  this  section,  and 
suppress  v]  the  statements  of  A1-A7  using  this  notation  should  be  clear. 
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•  Since  uji  =  iji  for  i  G  Z,  assumption  A1  follows  from  (4.8).  We  also  see 
that  the  sets  =  Si  and  =  Qi,  introduced  in  assumptions  A2  and 
A4  are  same.  Thus  A2  follows  from  (4.9). 

•  Assumption  A3  is  immediate  from  the  definition  ojt. 

•  Since  ojj  =  fjj,  the  set  introduced  in  assumption  A4,  is  given  by 

=  fli  =  Lix^eQifjj.  Since  r/j's  are  balls  of  radius  Rhj,  it  is  easily  seen, 
using  (4.18),  that  assumption  A4  is  satisfied  with  p  =  ZRC^ICy. 

•  RKP  shape  functions,  (f>j,  considered  here,  are  reproducing  of  order  fc  =  1, 
i.e.,  they  satisfy  (4.11)  with  k  =  1.  Thus  A5  is  satisfied  with  =  Ai  =  J 
(identity),  for  all  i  G  Z,  with  k  =  1. 

•  It  was  shown  in  [48]  that  if  the  weight  function  w(x)  €  C®,  then 

\\4>i\\w‘’°={f)i)  <  C{hi)~%  for  0  <  s  <  g,  for  i  G  Z. 

Thus  using  a  scaling  argument  and  this  estimate,  we  obtain 

\\<l>%\\H‘(f,i)  <  for  0  <  s  <  <7,  for  i  G  Z, 

where  h  and  hi  satisfy  (4.18).  Recall  that  we  assumed  <7=1. 

Xj  G  Qi-  Then 

and  combining  this  with  (4.21),  we  get,  for  0  <  s  <  1, 

ll'/’illffT’li)  —  C'(^i)~^'''”^^:  fo''  ^3  G  Qg  foi'  f  G  Z, 
which  is  assumption  A6  with  (7  =  1. 

•  A  scaling  argument  shows  that  assumption  A7  is  satisfied. 

We  remark  that  (3.61)  of  Remark  3.4,  together  with  condition  (b)  (following 

(4.18) ),  establishes  a  lower  bound  of  n,  namely  (n+  1)  <  k. 

We  have  thus  shown  that  assumptions  A1~A7,  with  fc  =  1  and  q  =  1,  are 
satisfied  by  RKP  particle-shape  function  systems  provided  (4.8),  (4.9),  (4.16), 

(4.18) -(4.20)  are  satisfied.  Thus  we  can  apply  Theorem  3.8  to  obtain  an  approx¬ 
imation  error  estimate  for  RKP  particle-shape  function  systems.  Note  that  the 
condition  (3.64)  in  Theorem  3.8  is  trivially  satisfied  with  A'^  =  /  for  all  x  G  X'' . 
We  remark  that  an  interpolation  error  estimate,  under  the  assumptions  (4.8), 
(4.9),  (4.16),  (4.18)-(4.20),  was  also  obtained  in  [48]. 

We  note  that  A1-A7  only  guarantee  good  approximability  of  the  shape  func¬ 
tions;  they  do  not  provide  a  recipe  to  construct  particle  shape  functions  that 
are  quasi-reproducing  or  reproducing  of  order  k.  In  fact  the  availability  of  such 
particle  shape  functions  is  assumed  in  A5.  Further  assumptions  may  be  needed 
to  construct  such  shape  functions;  for  example  (4.16),  (4.18)-(4.20)  were  needed 


(4.21) 
Now  let 
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to  construct  RKP  particle  shape  functions.  Therefore,  there  should  be  enough 
restrictions  on  the  particle  distributions  and  the  supports  of  shape  functions  so 
that  it  is  possible  to  construct  these  shape  functions  satisfying  A1-A7,  thereby 
ensuring  good  approximation  properties. 

Uniformly  distributed  particles: 

We  consider  the  uniformly  distributed  particles  Xj  =  jh,  j  G  Z”  as  in  Section 
3.2.  This  is  a  special  case  of  the  non-uniformly  distributed  particles  considered 
in  the  first  part  of  this  section.  For  each  Xj,  we  define  w^{x)  =  ),  where 

w(x)  >  0  is  a  continuous  weight  function  with  compact  support  r/  =  5^(0). 
Clearly,  r/j  =  supp  Wj{x)  =  B[ih{xj).  It  can  be  easily  shown  that  if  i?  =  3-y/n/2 
(in  fact,  we  need  only  R  >  y/n),  then  (4.8),  (4.9)  with  k  =  (4i?  +  1)”,  (4.18) 
with  Cl  =  C2  =  1,  and  (4.20)  with  C2  =  1/2  are  satisfied.  If  w{x)  =  w{r), 
with  r  =  ||a;||,  is  monotonically  decreasing  in  r,  then  it  also  can  be  easily  shown 
that  (4.19)  is  satisfied  with  Ci  =  w{^/n).  Therefore,  RKP  shape  functions 
(fiix),  associated  with  x/,  for  all  i  G  Z”  can  be  constructed  using  the  procedure 
described  in  (4.10),  (4.11)  and  (4.12)  for  k  =  1,  namely 

=  w^{x)  ^  (x  -  Xj^)“&(((x),  (4.22) 

|a|<fc 

where  {ba{x)}\a\<k  is  the  solution  of 

m((+^(x)6(((x)  =  5|^|,o,  \P\  <  k,  (4.23) 

\oi.\<k 

with  fc  =  1,  and 

m^(x)=  ^  u;^"(x)(x-x^")“.  (4.24) 

ieZ" 

Shape  functions  (pj’s  satisfy 

Y  P{Xj)(l>'}{x)  =  p{x),  for  all  x  G  M",  for  all  p  G  ■p'=(M”).  (4.25) 

As  with  the  non-uniformly  distributed  particles,  we  consider  the  family  of 
particle-shape  function  systems 

=  0<h<l, 

for  RKP  shape  functions  with  respect  to  uniformly  distributed  particles,  by  let¬ 
ting  =  {x  =  Xj  :  j  G  Z”}  and  using  h'f  =  h,  =  rj^  and  =  (f'f.  Note 
that  here  we  used  the  parameter  h  instead  of  v.  We  have  shown  above  that 
conditions  (4.8),  (4.9),  (4.18)-(4.20)  are  satisfied,  with  w{x)  =  w{r),  a  mono¬ 
tonically  decreasing  weight  function  in  r,  and  R  =  Therefore,  based 

on  the  discussion  on  RKP  particle-shape  function  systems  for  non-uniformly 
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distributed  particles,  it  clear  that  {A^^}o</i<i  satisfies  the  assumptions  A1-A7 
with  fc  =  1,  ensuring  good  approximation  properties  of  the  RKP  shape  func¬ 
tions. 

We  recall  that  in  Section  3.1,  the  particle  shape  function  was  defined 

in  (3.1)  by  scaling  and  translating  the  basic  shape  function  ^(cc)  for  uniformly 
distributed  particles,  i.e.,  they  were  translation  invariant.  We  will  show  that  the 
RKP  shape  functions  constructed  via  (4.22)  and  (4.23),  also  satisfy 

(3.1)  with  (f){x)  =  (j)Q{x)  {i.e.,  with  z  =  0  and  h  =  1).  We  assume  that  the  linear 
system  (4.23)  is  invertible  for  fc  >  1. 

From  (4.22)  and  (4.23)  with  z  =  0  and  h  =  1,  we  have 

(l){x)=w{x)  x°^bl{x),  (4.26) 

|a|<fc 

where  bl^{x)  are  the  solutions  of 

Y  =  ^|/3|.o,  for  \P\  <  k,  (4.27) 

and 

mi{x)  =  Yw{x-  j){x  -  j)“-  (4-28) 

jeZ" 

We  replace  x  by  in  (4.27)  and  (4.28)  to  get 


where 


,x  —  x^ ,,  ,x  —  x^ 


Y  ’^a+/3( — =  '^l/3fo>  for  1/^1  <  k,  (4.29) 


|q:|<A; 


^i(^)  = 


Yw{^^-j){^^-jy 


jeZ" 


jeZ' 

1 


h  h 


^  ^  w’jixyx-xY 


/zl“l 


((^)- 


Using  (4.30)  in  (4.29),  we  get 

El  fi  /  ~  \  r 

jT^m^+0{x)b^{——)  =  5|0|,o, 


|a|<A; 


and  therefore. 


Y  )  =  ^'^''^1/31.0  =  ^|/3|.o,  for  all  \(3\  <  k. 


|Q;|<fc 


(4.30) 


(4.31) 
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Since  &^(x)’s  are  unique  solutions  of  (4.23),  it  is  clear  from  (4.31)  that 


and  thus  from  (4.26),  we  have 


<t>{ 


X  — 


-)  =  w{ 


,X  —  Xj  -  ,  X 

h  ’ 

VI 

’  /il“l 

|a|<fc 

)“&o(a:) 

VI 

(j)’l{x). 

rh 

-) 


Thus,  for  uniformly  distributed  particles,  RKP  shape  functions  satisfy  (3.1), 
i.e.,  they  are  translation  invariant. 


Remark  4.1  To  approximate  functions  defined  on  a  bounded  domain  f2,  we  use 
the  restrictions  of  RKP  shape  functions  on  fl,  as  described  in  Section  3.3  (cf. 
(3.86)  and  Theorem  3.10).  We  note  that  the  RKP  shape  functions  corresponding 
to  the  particles  near  the  boundary  of  as  defined  here,  are  different  from  the 
RKP  shape  functions  defined  in  [48]  and  [58].  But  they  are  same  for  particles 
inside  sufficiently  away  from  the  boundary  dfl.  They  are  also  same  when 

n  =  M". 


4.2  Interpolation  and  Selection 

In  this  section,  we  will  address  the  interpolation  of  a  function  in  terms  of  particle 
shape  function,  and  will  propose  a  procedure  to  select  shape  function  that  will 
yield  efficient  approximation.  We  consider  uniformly  distributed  particles  {xj} 
in  M",  and  the  associated  particle  shape  functions  defined  in  (3.1),  where 
(j)  e  il^(K”)  with  g  >  1  has  compact  support;  supp  (j)  C  Br(0).  We  have  seen 
that  are  translation  invariant,  supp  (j)^  €  Bjih{Xj),  and  in  addition  they 
satisfy 

(4.32) 

We  assume  that  {(/'^}  are  reproducing  of  order  k,  i.e.,  (4.25)  holds. 

Let  n  be  a  bounded  domain  in  M".  We  will  consider  a  smooth  function  u(x) 
defined  in  and  study  the  error  u  —  IhU,  where  ThU  is  the  “interpolant”  of  u 
in  terms  of  (p’j.  The  results  in  this  subsection  are  from  [13],  and  we  refer  to  [13] 
for  some  of  the  details  that  we  do  not  present  here. 

We  now  define  the  “interpolant”  ThU  of  a  function  u.  For  any  x  G  M",  let 

A^  =  {kGZ-:xG 
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is  called  the  influence  set  for  the  point  x.  Then  {Xhu){x)  is  deflned  as 


{ihu){x)  =  ^  u{x^j)4)^{x).  (4.33) 


In  (4.33),  we,  of  course,  assume  that  u{Xj)  is  deflned  for  all  j  €  A^.  If  p  G 
then  from  (4.25)  we  have 

p{xj)(j)'^{x)  =  u{x^)(j)^{x)  =  p{x),  for  all  x  G  M",  (4.34) 

jeAJ  jez" 

i.e.,  IhP  =  P-  Now  let  u  G  H^{U)  with  s  >  nl2.  For  some  x  G  ft,  the  particles 
Xj  for  j  G  A'^  may  be  outside  fl,  and  u{xj)  may  not  be  deflned.  To  deflne 
Thu{x)  for  u  G  i/®(n)  and  for  all  x  G  fl,  we  need  an  extension  u  of  m  in  a  ball 
Bfig  containing  ft  such  that  dist{dQ,dBRg)  >  ph,  and  u  G  H^{Bjig).  Then, 

{ihu){x)  =  {ihu){x)  =  u{xj)(j)j{x),  for  X  e  (4.35) 

is  well  deflned.  For  an  extension  u,  we  may  use  u  =  Eu,  where  Eu  was  deflned 
in  (3.51).  Thus,  {Ihu){x)  for  x  G  fl  will  depend  on  few  values  of  u{xj),  where 
the  particle  Xj  is  just  outside  ft.  We  remark  that  IhU  is  not  an  interpolant  of  u 
in  the  usual  sense,  since,  generally,  (pjix'l)  yf  Sij,  and  hence  {Thu){xj)  yf  u{xj). 
We  deflne  the  function 

^a(x)  =  x“  —  {x'l)°‘ (pi  (x) ,  \a\  =  k  +  1,  for  x  G  M".  (4.36) 

i&A^ 

We  will  also  use 

^Q,(x)  =  Cq(x)  =  x“  —  i°‘(j){x  —  i),  \a\  =  k  +  1,  for  x  G  K”,  (4.37) 


where  A^  is  A^  with  h  =  1.  These  functions  will  play  an  important  role  in  the 
analysis  presented  in  this  subsection,  as  well  as  in  Section  5.  We  note  that  ^q(x) 
is  the  error  between  the  polynomial  x“,  with  |a|  =  k  +  1,  and  its  interpolant 
when  h  =  1.  In  1-d,  we  will  write  these  functions  as  ^^^^(x)  and  ^fe+i(x) 
respectively. 

We  begin  with  certain  results  about  these  functions.  We  first  present  some 
notations  that  will  be  used  in  these  results.  Let  /j*  be  the  cell  centered  at  Xj, 
deflned  by 

/(*  =  {x  :  ||x  -  x'^lloo  =  max  \xi  -  Xj,\ <  h/2}. 

For  each  Jj*,  we  deflne 


A';  =  {fc  G  Z"  :  77(1  n  y^  0}, 
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and 


We  note  that  cardinality  of  Aj  is  finite,  and  is  bounded  independent  of  j  and  h. 
Also  there  exists  ^  >  0,  independent  of  j  and  h,  such  that  Bj  C  B^  =  Bf^f^{x^) 
and  Li =  M”. 

Lemma  4.1  ([13])  ^a(x),  with  \a\  =  k  +  1,  is  periodic,  i.e., 

(^(x  +  x^)  =  e^(x),  for  any  x^].  (4.38) 

Proof.  We  first  note  that 

{x  +  Xj)"^  =  x‘^  +  p{x-,  x’)) ,  (4.39) 

where  p{x;  Xj)  is  a  polynomial  in  x  of  degree  <  k  with  coefficients  that  depend 
on  Xj.  Now  using  (4.39),  with  x  =  x'],  and  the  fact  that  the  (j)['s  are  translation 
invariant  and  reproducing  of  order  k,  we  get 

ieZ" 

ieZ" 

ieZ" 

=  E  +  E 

ieZ"  ieZ" 

=  (4.40) 

ieZ" 

From  from  (4.36),  (4.39)  and  (4.40),  we  get 

Cix+x'))  =  x--j2{x[rf>>[ix)  =  C{x), 

which  is  the  desired  result.  □ 

Lemma  4.2  ([13])  Let  a  =  a{i),  i  =  1,  •  •  •  ,  he  an  enumeration  of  the 
multi-indices  a  with  |q;(i)|  =  fc  +  1.  Let  Lj  be  the  cell  centered  at  the  particle 
Xj.  Then,  for  da  G  K,  we  have 

II  E  ^,daea{x)\\l,^,u)=h-^'^+-N^{A  +  K^B)N,  (4.41) 

|a|  — fc+1 


ieZ" 


43 


where  V  =  ^0(2),  •  ■  ■  >  da{Mk)Y  A,  B  are  Mk  x  matrices  given  by 

1 

a{l)la{m)l 
1 

a{l)la{m)l' 


Aim  — 
Blm.  — 


I  '  ^^a(m)  dXj 

rCa(/)^a(m)  dx  , 


(4.42) 

(4.43) 


respectively,  and  I  =  [—1/2, 1/2]". 


Note:  The  matrices  A  and  B  are  independent  of  Ij. 

Proof.  A  simple  scaling  argument,  used  with  (3.1),  shows  that 

Now,  using  the  periodicity  of  ^a(x),  a  standard  scaling  argument,  and  this 
identity,  we  have 


|2 


|a|=fc+l 


|2 


|Q;|=fe  +  l 


=  h2('=+i)h"-^|j  ^  -,d^yUy)\\HO(i) 


|a:|=fe+l 


=  h2'=+"V^AV. 


(4.44) 


Using  a  similar  argument,  we  have 

1 


\a.\—k-\-l 

Combining  this  identity  with  (4.44),  we  get 

|Q:|^fc+l 

which  is  the  desired  result.  □ 


Lemma  4.3  ([13])  Let  ij*  he  the  cell  centered  at  the  particle  x^ ,  and  consider 
the  corresponding  set  Bj.  Suppose  u  G  i7^+^+^(i3j*)  with  <?  >  f  when  n  >  2, 
and  q  =  0  when  n  =  1 .  Then, 

(a)  for  any  (5  >  0, 

<  (1  +  52)11  E  ^P“«)(^7)eo(^)ll?7i(7.) 

+(1  +  ^)Ch2'=+2  ^  \\D^u\\l^^^,y  (4.45) 

|Q:|=fe+2  ^ 
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and 

(b)  for  any  5  >  0, 

II  <  {l  +  S^)\\u-i^u\\l.^,,^ 

|a|— fc+1 

+  {^+ \\D°‘u\\lj^^g^y 
|a|— fc+2 

(4.46) 

The  proof  of  this  result  is  based  on  Taylor’s  Theorem,  a  bound  on  the  remain¬ 
der  in  Taylor’s  Theorem,  and  a  bound  on  the  interpolant  of  the  same  remainder. 
We  do  not  include  the  proof  here,  and  refer  to  [13]. 

We  will  now  study  the  interpolation  error  u  —  ThU,  where  u  is  a  smooth  func¬ 
tion  in  f2.  An  interpolation  error  estimate,  namely  ||m— 2'^M||_f/i(a)  «  0{h^),  was 
proved  in  [48,  59]  for  the  RKP  shape  functions.  A  similar  order  of  convergence 
in  the  norm  was  also  obtained  for  MLS  shape  functions  in  [1,  2].  We 

note  that  the  definitions  of  XhU  for  the  RKP  shape  functions  and  MLS  shape 
functions,  presented  in  these  papers,  are  slightly  different  from  our  definition 
as  given  in  (4.35).  From  the  proof  of  the  our  next  result,  we  will  obtain  an 
estimate  of  ||u  —  Ihu\\H^{n)  where  the  shape  functions  shape  functions  are  re¬ 
producing  of  order  k.  Moreover,  this  theorem  gives  some  information  on  the 

size  of  which  facilitates  the  selection  of  “good”  shape  functions, 

which  will  be  discussed  later. 

We  now  present  the  main  result  of  this  section.  We  define  certain  sets,  which 
will  be  used  in  this  result: 

A'*  =  {A:  G  Z"  :  L!  n  0}, 

=  {A:  G  Z"  :  4^  C  fl},  0^  = 

It  is  clear  that  C  14  C  and  [H  —  fi^|  ^  0,  [flft,  —  14|  — >  0  as  A-  — >  0.  Also 
nc  C  and  \lt  -  14|  ^  0,  |i?^  -  14|  ^  0  as  Ai  ^  0. 

Theorem  4.1  ([13])  Let  A  he  the  largest  eigenvalue  of  the  matrix  A  given  in 
(4.4^)'  Suppose  9  >  f  when  n>2,  and  q  =  0  when  n  =  1.  Then,  we  have 

||u  — 

sup  lim  -  .  - =  A, 

„g//fc+2+<,(o)  ^^0  h^'^Qh{u) 

where 

Qh{u)  =  \u([jk+l(^Q^-^+  h  'Y  l|-D“M||ffg(0)- 

|a|— fc+2 


(4.47) 

(4.48) 
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Note:  In  (4.47),  we  consider  u  G  such  that  u  ^  7^*. 

Proof.  We  will  first  prove  that  for  u  € 


lim 

h^O 


||w 

K^’^Qhiu) 


V'^{x)AV  (x)  dx 


(4.49) 


where 

V'^{x)  =  [D°‘^^'^u{x),D°‘^'^\{x), ...,  D“(“'‘)u(x)], 

and  a{i),  1  <  i  <  Mk,  are  the  multi-indices  with  |Q;(t)|  =  k  +  1. 

Let  u  G  77^+^+"? (n),  and  suppose  u  is  an  extension  of  u,  as  discussed  before. 
Since,  fl  C  fl/i,  we  have 


llu  —  11"*^  “  'y  11^  ^hu^fji(^ihy 

j&A>' 


Therefore,  using  (4.45),  (4.41),  and  recalling  that  =  Dj^^hBj,  we  get  for 
any  d  >  0, 

<  (1  +  <5^)E|I  E  ^(^“^)(^")5a(^)ll^i(Z.) 

|a|=fc+l 

+(i + i)ck“+>  y:  y: 

JGA’'  |a|  =  fc-|-2  ■’ 

<  {l  +  6‘^)h‘^'^  h^V[{A  +  h‘^B)Vj 

jeA’' 

+(1  +  1)C72'=+2  ^  (4-50) 

|a|=fe-|-2 


where 

Therefore,  dividing  (4.50)  by  hf^Qh{u),  where  Qh{u)  is  defined  in  (4.48),  we  get 


h‘^^Qh{u) 


Qh{,A) 


+  + 


1  ,  2  5Z|Q|=fe+2  11-^ 


Qh{A) 


(4.51) 


A  typical  term  of  the  quadratic  form  {A  +  h?B)Vj  is 
D^^'^u{x']){Aa  +  h^Bu)D^^^h{x'^). 

Since 

Im  ^  /i”L'“WM(a;j^)A,,T»“(')u(x^^)  =  f  D°^^'^u{x)AiiD°^^^^u{x)dx 


jeA'' 
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and 


we  have 


lim  y  h'^D°‘^^u{x'])BaD°‘^^'^u{x'])  =  0, 
h^O  ^  ^  ^ 


jeA^ 


lim  y  h'^V^{A  +  h‘^B)Vj=  f  V'^{x)AV{x)dx. 


jGA^ 


(4.52) 


Since  \B^  —  n|  ^  0  as  ft-  ^  0,  we  have 

lim  y  =  J2  (4-53) 

|a|  — fc+2  |q:|— fc+2 

Also  liuih^oQhiu)  =  |M|_f/<=+i(o).  Thus,  for  any  (5  >  0,  using  (4.52)  and  (4.53) 
in  (4.51),  we  get 


lim  sup 


|u- ^  ^  2AnV^(x)AV(x)dx 


h^o  d^'^Qhiu) 
and,  since  ft  >  0  is  arbitrary,  we  have 


<  (l  +  ft")  = 


J^V^{x)AV{x)  dx 

Jr  h^'^Qf.iu)  -  ■ 


(4.54) 


Following  the  argument  leading  to  (4.54),  but  using  A^,  B^,  and  (4.46)  instead 
of  A^,  B^,  and  (4.45),  respectively,  we  can  also  show  that 


V'^{x)AV (a;)  dx 


<  lim  inf 


\u-ihu\\jji 


(n) 


(4.55) 


i^i//fc+i(Q)  h—*o  h‘^^Qfi(u} 

Combining  (4.54)  and  (4.55),  we  see  that  lim^^o  exists,  and 

ll'w- J'wllffi(O)  _  f^V'^(x)AV(x)dx 


tz  h^kQAu) 


which  is  (4.49). 

Since  A  is  the  largest  eigenvalue  of  the  matrix  A,  from  the  usual  variational 
characterization  of  eigenvalues,  we  have 


Mfc 


f  V^{x)AV{x)  dx  <  X  [  '^\D°‘^^'^u{x)\'^  dx  =  X\u\\jk+i/Q;y 
Jo.  Jn  ^ 

Thus  from  (4.49)  we  get 

tormy„eH‘+^+«(n). 
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Hence 


sup  lim 

ueHP+'^+i{Q) 


llw 

h'^^Qhiu) 


<  A. 


(4.56) 


Let  V  =  [vi,V2,‘"  be  an  eigenvector  of  A  corresponding  to  A.  Then  it 

is  easily  seen  that  there  is  a  m  G  7^^+^  such  that  the  vector  V(x)  =  v.  For  this 
particular  u,  we  have 

J^V^{x)AV{x)  dx  ^  - 


Hence,  from  (4.56)  we  conclude  that 


sup  lim 

„g_f/fc+2+5(0)  h^o 


llu  ^huW'n'^in) 
K^’^Qhiu) 


=  A, 


which  is  the  desired  result.  □ 

Remark  4.2  We  know  from  (4.35)  that  the  interpolant  of  a  smooth  function 
depends  on  its  extension  to  M”.  But  it  is  clear  from  the  proof  of  Theorem  4.1 
that  (4.47)  is  valid  for  any  extension  satisfying  (3.51). 

Remark  4.3  We  note  that  same  result  holds  for  the  i7^-seminorm  of  the  in¬ 
terpolation  error,  i.e.,  for  q  >  ^  when  n  >  2,  and  q  =  0  when  n  =  1,  we 
have 

sup  lim  „ - — = - — - — ^ - r  =  A 

^^0  h  +  h22i\a\=k+2 

Remark  4.4  From  (4.51)  in  the  proof  of  Theorem  4.1,  we  can  obtain  an  inter¬ 
polation  error  estimate. 


\\u  —  ihu\\m(Q:)  <  C'^^l|w||//fc+2+<,(n), 


where  C  may  depend  on  12,  but  is  independent  of  u  and  h.  We  note  however, 
that  this  is  not  the  optimal  error  estimate.  For  an  outline  of  the  proof,  see  [13]. 

We  have  seen  in  Remark  4.4  that  if  the  particle  shape  functions  are  repro¬ 
ducing  of  order  k,  then  for  a  smooth  function  u, 

\\u-Xhu\\mm-o{h’^), 

where  ThU  is  the  interpolation  of  u  as  defined  in  (4.35).  There  are  many  classes 
of  shape  functions  that  have  these  properties.  We  have  seen  in  Section  4.1 
that  translation  invariant  RKP  shape  functions  depend  on  the  weight  function 
w{x),  and  different  choices  of  w{x)  will  generate  different  classes  of  such  shape 
functions. 
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We  will  assess  the  aproximability  of  a  family  }  of  shape  functions  by 
the  size  of  A,  the  largest  eigenvalue  of  the  matrix  A  defined  in  (4.42).  We  note 
that  A  is  computable,  and  depends  only  on  the  basic  shape  function  (j){x).  We 
emphasize  that  A  does  not  depend  on  u  or  on  h.  From  (4.47),  we  know  that 


\\u-ihu\\Hi{n)  < 
hP^/Qhiu)  "" 


for  small  h. 


Thus  we  see  that  A  is  a  useful  measure  of  the  approximabilty  of  the  family  {(/>j  }, 
determined  from  the  basic  shape  function  ^(x). 

We  will  illustrate  our  selection  scheme  in  1-d,  and  will  rank  the  shape 
functions  according  to  to  their  approximability.  We  note  that  in  1-d,  A  = 

( ^^'“(^+1)!°’^*  )^-  the  rest  of  this  paper,  we  will  suppress  77^(0, 1)  in  |C/c+i|ffi(o.i)> 
and  instead  write  |^fe+i|i. 

We  considered  four  different  classes  of  RKP  shape  functions,  reproducing  of 
order  1,  corresponding  to  four  different  weight  functions  w{x).  These  w(a;)’s 
were  (4.2)  with  6  =  2,  (4.3),  and  (4.4)  with  I  =  2,4.  We  then  computed  |^fc+i|i 
for  each  of  these  four  classes  of  shape  functions  for  R  =  1.7;  we  obtained 


0.237,  for  w{x)  in  (4.4),  I  =  2 
0.203,  for  w{x)  in  (4.2),  5  =  2 
0.095,  for  w\x)  in  (4.3) 

0.029,  for  w\x)  in  {4.4), I  =  4 


We  choose  the  RKP  shape  functions  corresponding  to  w{x)  given  in  (4.4)  with 
I  =  4,  since  these  shape  functions  yield  the  smallest  value  of  |^fe+i|i.  We  note 
that  the  value  of  |Cfe+i|i  depends  strongly  on  R,  and  the  shape  function  cor¬ 
responding  to  w{x)  given  in  (4.4)  with  I  =  4  may  not  be  our  choice  for  other 
values  of  R.  We  refer  to  [13]  for  futher  discussion  on  this  issue. 

To  validate  our  criterion  of  selection  of  the  shape  functions,  we  have  con¬ 
sidered  the  function  u{x)  =  x'^  on  the  interval  ft  =  (0, 1)  and  computed  the 
error  |m  —  ihU  is  the  interpolant  of  u  with  respect  to  the  four 

classes  of  RKP  shape  functions  described  in  the  last  paragraph,  with  h  =  1/n, 
n  =  40,  50, ... ,  100.  We  note  that  the  definition  of  IhU  requires  the  values  of 
u{x)  in  a  small  neighborhood  of  fl,  and  we  have  extended  u  =  x^  outside  12  by 
itself.  We  summarize  the  results  in  the  following  table. 
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n 

\u  —  ZhU\H^(Q) 

Conical:  1  =  2 

Gauss:  <5  =  2 

Cubic  Spline 

Conical:  1  =  4 

40 

1.607e-2 

1.376e-2 

6.435e-3 

2.283e-3 

50 

1.281e-2 

1.096e-2 

5.130e-3 

1.730e-3 

60 

1.066e-2 

9.112e-3 

4.267e-3 

1.396e-3 

70 

9.126e-3 

7.800e-3 

3.653e-3 

1.172e-3 

80 

7.980e-3 

6.819e-3 

3.194e-3 

1.012e-3 

90 

7.090e-3 

6.058e-3 

2.838e-3 

8.908e-4 

100 

6.379e-3 

5.449e-3 

2.553e-3 

7.962e-4 

Table  4.1:  The  iif^-seminorm  of  the  error,  \u  —  where  ThU  is  the 

interpolant  of  u{x)  =  a:^  using  RKP  shape  functions  that  are  reproducing 
of  order  1,  corresponding  to  different  weight  functions  w(x).  The  radius  of 
support  of  oj(x)  is  R  —  1.7. 

From  Table  4.1,  it  is  clear  that  the  error  |u  —  ThM|_f/i(n)  can  be  ranked 
according  to  the  size  of  |^2|i  for  the  four  choices  of  w(x)  considered  here  with 
R  =  1.7;  the  error  and  |^2|i  are  both  minimal  when  w{x)  is  the  conical  weight 
function  with  I  =  4. 

This  selection  scheme  is  based  on  (4.47),  and  we  know  from  Remark  4.2  that 
(4.47)  is  valid  for  any  extension.  We  refer  to  [13]  for  an  experimental  illustration 
of  this  fact.  We  remark  that  this  selection  scheme  is  also  valid  for  the  projection 
error,  which  will  be  indicated  by  our  results  in  the  next  section. 

5  Super  convergence  of  the  gradient  of  the  solu¬ 
tion  in  L2 

Superconvergence  is  an  important  feature  of  finite  element  methods,  which  al¬ 
lows  an  accurate  approximation  of  the  derivatives  of  the  solution  of  the  under¬ 
lying  BVP.  In  this  section,  we  will  discuss  the  idea  of  superconvergence  when 
particle  shape  functions  are  used  to  approximate  the  solution  of  a  BVP.  We  will 
consider  uniformly  distributed  particles  and  the  associated  particle  shape  func¬ 
tions,  which  were  developed  in  Sections  3.1  and  3.2.  For  uniformly  distributed 
particles,  a  careful  analysis  in  1-d  can  be  easily  generalized  to  higher  dimen¬ 
sions.  Thus,  in  this  section,  we  present  the  results  in  1-d,  but  by  restricting  our 
analysis  to  1-d,  we  avoid  some  details  that  arise  in  higher  dimensional  analysis. 

We  will  use  the  notation  that  was  introduced  in  Section  3.1,  but  restricted 
to  1-d,  i.e.,  for  ft,  >  0,  we  consider  x’j  =  jh,  j  G  Z,  and  the  corresponding 
shape  function  defined  in  (3.1).  We  assume  that  the  shape  functions  are 
reproducing  of  order  k.  We  use  the  following  notation: 

=  {x^,  A';  =  {meZ:r,tnl^^  0}; 

^i  =  (i.J  +  l),  (with  ft  =1). 
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We  assume  that 
or  equivalently, 


card(Aj)  <  k, 
card(A^)  <  K, 


where  k  is  independent  of  j  and  h.  We  assume  that  the  basic  shape  function 
(j){x)  is  such  that,  for  any  v{x)  =  ^  ^  ^0,  there  exist  positive 

constants  Ci,  C2,  independent  of  v,  but  may  depend  on  k,  such  that 


Cl  c^^  <  j  dx  <  C2  ^  cl  (5.1) 

jeAo  jeAo 


This  implies  that  the  functions  {(j)i{x)}i^Ao  Etre  linearly  independent  in  /q,  i.e., 


Cj(j)j{x)  =  0,  X  G  Iq  implies  cj  =  0,j  G  A^. 

j&Ao 


Throughout  this  section,  we  use  C,  Ci,  C2  as  generic  constants,  which  will  have 
different  values  in  different  places. 

Consider  =  (— c,  d)  C  K.  Let  uq  G  be  the  solution  of  the  Neumann 

problem 

B(uo,v)  =  T(v),  for  all  u  €  (5.2) 

where 

B{u,v)  =  /  {u'v'  +  uv)dx  and  B{v)  =  /  fvdx 

Jn  Jn 

as  in  (2.4)  and  (2.5).  We  will  often  use  the  notation  B^{u,v)  to  denote  the 


above  bilinear  form,  where  the  12  is  replaced  by  another  domain  F. 

Let  Uh  G  be  the  solution  of 

B{uh,v)  =  F{v),  for  all  V  G  (5.3) 

where  was  defined  in  (3.55).  It  is  clear  from(5.2)  and  (5.3)  that 

B{uo  -  Uh,  v)  =  0,  for  all  v  G  (5.4) 

and  we  easily  have 

l|w?i||ii-i(n)  <  ||uo||hi(o).  (5.5) 

Recall  that  the  functions  in  are  restrictions  of  the  functions  in  Sh  = 
on  12  (c/.  (3.2)  and  (3.55)).  Thus  (5.4)  is  true  when  is  replaced  by  Sh- 
We  assume  that  for  any  p  >  0, 

11^0  —  UhllL2(Bp(0))  <  C/l^''‘^||'Uo||_f/fc+i(0)p2  ,  (5.6) 


and  there  are  positive  constants  Ci,C2,  independent  of  u,  h,  and  p,  such  that 


Cih'=p5  < 


ko  -<IU2(B,(0)) 


B'=+i(n) 


<  C2/i>C 


(5.7) 
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where  Bp{0)  =  {x  :  \x\  <  p}  and  Bp{0)  C  fl.  We  will  write  Bp  =  Bp{0)  through 
out  this  section. 

The  main  goal  of  this  section  is  to  investigate  the  error  u'{x)  —  u'f^{x)  in  a 
neighborhood  of  x  =  0,  i.e.,  for  x  €  Bfj  CC  and  H  =  W ,  7  <  1,  where  7  will 
be  chosen  later.  We  will  prove  the  following  result: 

Theorem  5.1  Suppose  uq  and  Uh  satisfy  (5.6)  and  (5.7),  and  let  Ch  =  uq  —  Uh. 
Moreover,  assume  that  uq  G  W^^{B2h)  ■  Then  for  h  small  enough,  there  exists 
e*  >  0,  such  that 

\Wh  -  r('«o)?fc+/||L2(iJ„)  ^  e* 

\K\\lm 

where  T{uo)  =  and  Cfc+i(a^)  =  fk+i  is  defined  in  (4.37). 

Remark  5.1  Theorem  5.1  is  a  superconvergence  result.  It  shows  that 
114  -  ^('Wo)Cfe+l  lk2(Bff)  <<  I|4IIl2(Bh)- 

This  allows  one,  for  example,  to  analyze  the  effectiveness  of  an  error  estimator 
as  was  done  in  [19]. 

Since  the  all  the  results  in  this  paper  have  been  presented  in  terms  of  L2 
based  norms  {i.e.,  in  terms  of  the  usual  Sobolev  norms),  we  also  present  this 
result  in  terms  of  L2  based  norm.  Superconvergence  in  L^o  will  be  addressed  in 
a  forthcoming  paper.  Assuming  superconvergence  in  Loo,  the  superconvergence 
points  and  superconvergence  recoveries  in  the  case  of  particle  shape  functions 
can  be  obtained  analogously  as  in  [19].  At  the  end  of  this  section,  we  will  see 
an  example  where  the  superconvergence  points  are  distributed  differently  than 
in  the  classical  FEM. 

Remark  5.2  The  essential  aspects  of  superconvergence  analysis  in  the  classical 
FEM  are  interior  estimates,  developed  in  [71],  [75],  [88].  This  analysis  strongly 
utilizes  the  polynomial  character  of  the  shape  functions.  Here,  in  the  case  of 
particle  shape  functions,  we  had  to  develop  another  approach  to  the  analysis 
of  superconvergence,  which  is  based  on  weighted  Sololev  spaces.  The  main 
idea  of  the  proof  of  our  superconvergence  result  is  to  show  that  locally,  the 
approximation  error  is  asymptotically  same  as  the  error  in  the  interpolation 
of  a  polynomial  of  degree  fc  +  1  by  particle  shape  functions.  The  analysis  is 
technical;  we  present  the  main  idea  of  this  analysis  in  this  section. 

Remark  5.3  Assumptions  (5.6)  and  (5.7)  are  directly  related  to  the  control  of 
pollution,  as  in  FEM.  The  assumption  that  uq  G  H4>+^(Lp)  is  analogous  to  the 
assumption  in  FEM  (see  [19]). 
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To  prove  Theorem  5.1,  we  will  first  develop  certain  ideas  and  establish  many 
technical  results.  Towards  this  end,  for  given  parameters  H  =  ,  with  7  <  1, 

and  a  >  1,  we  define  the  function  g{x)  by 

1,  -H  <x<H 

x>H  (5.8) 

^a{H+x)^  x^-H. 

where  a  is  such  that  ah  <  1,  and  will  be  chosen  later.  We  note  that  a  proper 
choice  of  7  and  a  is  crucial  for  the  analysis  presented  in  this  section.  Often,  we 
will  use  g  =  g{x),  gi  =  g{x’^)  and  g.+  i  =  g{x’^  +  h/2). 

Generalized  Interpolant  and  certain  norm  estimates: 

We  first  introduce  the  idea  of  generalized  interpolant  of  a  function  u,  which 
is  different  than  defined  in  Section  4.2.  Let  Iq  =  I-i  U  /q  U  {0}  =  (—1, 1) 
and  Aq  =  A-i^JAq.  Then  from  (5.1),  it  is  clear  that  there  are  positive  constants 
Cl,  C2,  independent  of  u  =  '^i^zCi(j)i{x),  but  may  depend  on  k,  such  that 

Cl  cl  <  f  dx  <  C2  ^  c)',  (5.9) 

iedo  jGAo 


which  implies  that  {(f>i{x)}^^^^  are  also  linearly  independent  in  /q.  We  define 
ipoix)  =  0'i4>i{x)  with  supp  'ipo  =  lo  (closure  of  Iq),  such  that 


Ho 


^o{x)(l)o{x)  dx 
'ipo{x)(l)j{x)  dx 


Using  (5.9),  we  can  show  that 


1, 

0,  for  all  j  G  Ao,  j  A  0- 


(5.10) 


Wi^ohoHo)  <  c- 


(5.11) 


We  also  note  that,  since  form  a  partition  unity  on  Iq,  from  (5.10) 

we  have 


Co 


Ho 


E 

Ao 


ijj^{x)dx=  /  ■0o(^)  =  1. 


(5.12) 


Ho 


Let  Aiix)  =  V'o(f  ~  *)•  Then  supp  Ai  = 
Note  that  =  M.  Now  for  a  given  v  G  L2(®)- 

interpolant  of  v  as 

ihv{x)  =  ^4'f(u)(()f(a:), 

iCil 


where  =  {xl_^,xY). 
we  define  the  generalized 

(5.13) 


53 


where 


(5.14) 


^  'ipi{x)v{x)  dx. 

We  note  that  T^v{x)  depends  on  the  v{y)  for  y  G  where  A^{x)  = 

G  Z  :  x  G  rjj^}.  We  also  define 

Lemma  5.1  Suppose  v(x)  =  X)i6Z 

=  ^^(v)  (5.15) 

and  I^v(x)  =  v(x).  (5.16) 

Proof.  From  (5.10)  and  the  definition  of  '!'(*(«)  in  (5.14),  i  G  Z,  we  have 
4'?(w)  =  ^  J^^f;'f{x)v{x)dx 

•  jed'* 


which  is  (5.15).  Now  using  (5.15)  in  (5.13).  we  get  (5.16).n 

Remark  5.4  We  note  that  if  ri  is  a  linear  combination  of  (f^’s  only  locally, 
i.e.,  in  a  bounded  open  interval,  then  =  v  only  in  the  interior  of  that  open 
interval.  More  precisely,  T^v  =  u  in  an  interval  /  if  u  is  a  linear  combination  of 
(fi  S  in  Ux^I  UigA'*(a;)  ■ 

We  will  use  the  following  result  later. 

Lemma  5.2  Let  LI  he  a  hounded  interval,  and  suppose  u  G  L2{Ll).  Then 

\\ilEu\\m(K)  <  Ch~^\\u\\L^(Q,) 

where  E  is  the  extension  operator  satisfying  (3.51). 

Proof.  We  first  note  that  the  extension  Eu  of  u  satifies  \\Eu\\l.^i^^')  < 
C'||M||L2(a)-  Now,  from  (5.13)  and  (4.32) 

jeA) 

<  Ch-^  Y.  (5-17) 

jeA^ 
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where  C  depends  on  k;  and  using  Schwartz  inequality  on  (5.14)  with  v  =  Eu, 
and  a  scaling  argument,  we  get 

^3  ^3 

<  l\ml^ijJ^^{Eu)^dx].  (5.18) 

^3 

Thus,  from  (5.17),  (5.18),  and  the  fact  that  ||ifM||i2(R)  — 

W^hEuW'jji^g^  <  Ch  ^||'u||i2(n)) 
which  is  the  desired  result.  □ 

Remark  5.5  We  can  also  show  that 

\\ihEu\\L2iW)  <  <^11^1^2(0) 

using  the  same  arguments  as  in  the  proof  of  Lemma  5.2. 

Consider  the  function  v{x)  =  X^ieA'*  o^i  Then,  using  scaling, 

translation,  and  (5.1),  we  have 

Cih  {c'lf  <  f  v‘^dx<  C2h  Y  (5.19) 

where  Ci,  C2  are  positive  constants,  independent  of  h  and  j,  but  may  depend 
on  K.  Using  (5.19),  we  can  show  that  if  v{x)  =  =  0  in  L2,  then 

=  0,  for  all  z  G  Z,  i.e.,  {(j)^}  are  linearly  independent. 

We  will  now  prove  certain  lower  bounds  for  gv^  dx  and  gv'  dx,  where 

^3  ^3 

g{x)  has  been  defined  before.  We  first  prove  the  following  inequality. 

Lemma  5.3  Let  zq,  zi  he  integers  such  that  zq  <  zi,  and  suppose  {ci}lE^  are 
real  numbers.  Then  there  exists  a  positive  constant  C,  depending  only  on  zi  —  zq, 
such  that  for  any  k,  io  <  k  <  ii,  we  have 

il  il-1 

^5i+i(c*  -Cfc)^  <EY  g^+^ici+l  -  c,f.  (5.20) 

i—io  i=io 

Proof.  Suppose  the  integers  zo,zi  are  such  that  H  <  i^h  <  iih,  where  H  = 
h'* ,  7  <  1.  Then 

il  k—1  il 

^5i+i(c* -Cfc)^  =  ^g,+  i(ci  -  Cfe)2  +  Y  9z+i{ci- Ck)^-  (5.21) 

i—io  i—io  i—k-\-l 
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We  first  note  that 


9,+  iici-Ckf 

ii  i—1 

i—k-\-l  j—k 

i—k-\-l  j—k  \ 


9j+h 


But  from  the  definition  of  g{x)  in  (5.8),  we  have 


(1 


9j+i 


and  using  this  in  (5.22),  we  get 

ii  ii  i—l 

Y  9^+^{ci-Ckf  <  C  Y  H5i+i(cj+i -Cj)" 

i—k-\-l  i—k-\-l  j—k 

ii-1 

—  ^  5j+i  (Cj  +  l  ~  Cj)  , 

j=k 

where  C  depends  on  {ii  —  ig). 

Using  similar  arguments  we  can  show  that 


fe-i 


fe-i 


(5.22) 


(5.23) 


(5.24) 


Y9^+k^'^k- Cif  <  C{k-l-io)Y9j+k{cj+i- Cjf,  (5.25) 


where  C  depends  on  (U  —  ig).  Therefore,  combining  (5.21),  (5.24),  and  (5.25) 
we  have 

n  ii 

-Cfc)^  <CY9j+i{c3+i  -  Cjf,  (5.26) 

i—iQ  i—io 

where  C  depends  on  (z*  —  zq).  Using  similar  arguments,  we  can  prove  (5.26)  for 
all  integers  zg,ii  such  that  ig  <  zi.  □ 


Lemma  5.4  Suppose  v{x  )  =  Then 

(a)  there  are  positive  constants  Ci,C2,  independent  of  v,  h  and  j,  but  may 
depend  on  k,  such  that 

Cih  Y  9^{c^y  <  [  gv"^  dx  <  C2h  Y  9t{c^fi  (5-27) 

ieA'*  ^3 

3  3 
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(b)  there  is  a  positive  constant  C,  independent  of  v  and  h,  such  that 

\  5i+i(c*>i  -  <C  [  gv'^dx.  (5.28) 

Proof,  (a)  Consider  j  G  Z  and  the  corresponding  A’f  such  that,  for  i  G  Aj, 
H  <  x^.  Let  gM  =  max^g^h  and  g^  =  min^g^h  g{x^).  Then,  it  is  easy  to 
check  that  ^  <  C,  where  C  depends  k.  Now,  using  (5.19),  we  have 


ieA^ 


< 


< 


gnh  ^  {c’lf 

ieA^ 


9m  / 
Cigm  J 


gv^  dx 


<  Cl  gv^  dx. 


(5.29) 


Using  a  similar  argument,  we  get 

[  gv^dx  <  C  9i{c^f 


ieA’' 


Combining  the  above  with  (5.29)  gives  the  required  result.  Using  similar  argu¬ 
ments,  we  can  prove  (5.27)  for  any  j  G  Z.  □ 

(b)  Let  u  =  '^i^j^Ci(j)i{x).  Then  from  (5.14)  and  (5.15)  with  /i  =  1,  we  have 


Ci  = 'I')('u)  =  /  'tpl{x)u{x)  dx, 
Ji-l 


(5.30) 


and  therefore. 


r^+2 


Q  +  l  Q  — 


'4^ij^\{x)u{x)  dx  —  /  'ipl{x)u{x)  dx 

Ji-l 


ni-\-2 


-  ipl{x))u{x)  dx. 


(5.31) 


Let  F{x)  =  J^_i['4’i+i{t)  —  'ilj}{t)]dt.  Using  translation  and  (5.12),  it  is  easily 
seen  that 


/  i’lit)dt=  V’i+i(i)rfi=l 
Ji-l  Ji 


and  therefore,  F{i  —  1)  =  F{i  -|-  2)  =  0.  Also,  using  the  Schwartz  inequality  and 
(5.11),  we  can  show  that 


^z+2 


F"  dx  <  C. 


Ji-l 
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Now,  using  the  above  bound,  integrating  (5.31)  by  parts,  and  using  the  Schwartz 
inequality,  we  get 


(ci+i  —  Ci)^  =  (  /  Fu'  dxY  ^  C'  /  u'^  dx. 

Ji-l  Ji-l 


(5.32) 


Let  V  =  Then  by  a  standard  scaling  argument,  we  have 


r^i+2  1 

{v'{x)fdx=-l  {u{y))'^dy, 

-rh  fl 


(5.33) 


li-l 


where  u{y)  =  Therefore,  from  (5.32)  and  (5.33)  we  have 

..h 

1 


^(cti-c'^)^<C 


dx. 


(5.34) 


From  the  definition  of  g{x),  we  can  show  that  for  x  €  {x^_i,x^j^2)j 
(1  +  ^  C'-  Therefore, 


^5,+  i(cf+i-cf)^  <  Cl  gv'\l  + 


Fn  :  9^+h-9 


<  C 


C^i  +  2 


9 


)  dx 


gv'^  dx, 


and  hence, 

^^9,+  ii4+i- <  C^J^  gv'"^  dx  <  C  J^gv'^  dx, 


which  is  the  required  result.  □ 


Remark  5.6  We  note  that  it  is  possible  to  show  that 

f  gv'^dx  <C^Y1  gz+i{ci+i  -  cf)^ 

and  together  with  (5.28)  we  see  that  ^  9i+^i'^i+i  ~  is  equivalent  to 

l^lffi(R)-  proof  of  this  fact  is  easier  than  the  proof  of  (5.28),  and  we  do  got 
provide  the  proof  here. 


A  perturbed  bilinear  form  Bq{u,v)  and  related  results: 
For  a  given  0  >  1,  we  now  consider  the  bilinear  form 

B^{u,  v)  =  B®(u,  v)  +  0£>“(m,  v). 
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where 

B^{u,v)  =  /  {u'v'  +  uv)dx  and  D^{u,v)  =  /  uvdx. 

Jw  Jw 

We  will  write  Bq{u,v)  =  Bq{u,v),  but  will  use  Bq{u,v)  when  the  domain  of 
integration  is  F  instead  of  M.  Also  we  will  use  D^{u,  v),  where  M  is  replaced  by 
a  domain  F  in  the  definition  of  D^{u,v). 

Let  Hg  Q  and  q  be  Hilbert  spaces  defined  as 

HgQ  =  {m:||'u||i  0=  /  gu^dx  +  {l  +  Q)  /  gu^  dx  <  oo}; 

Jr  Jr 

^g-^e  =  :  llwllpg-ye  =  /  g~^u^  dx  +  {1  +  Q)  (  g~^u^  dx  <  oo}. 

IR  M 

We  will  choose  0  later.  The  choice  of  0,  along  with  the  choices  of  7  and  a, 
mentioned  before,  is  important  for  the  main  result  of  this  section.  We  assume 
that  %-  <  1  where  0=1  +  0. 

We  will  often  suppress  0  in  ||M||i,g,e  and  HuHi^g-ye  and  instead  write  HuHi^g 
and  ||u||i_g-i  respectively.  We  will  also  use  the  fact  that  \g' /g\  <  a,  which  is 
obvious  from  the  definition  of  g{x). 

Remark  5.7  The  space  q  is  directed  towards  obtaining  interior  estimates 
of  e'^,  i.e.,  ejj  is  locally  characterized  through  the  use  of  the  space  Hg  Q. 

We  now  consider  Bq{-,  •)  :  q  x  q  ^  K. 

Lemma  5.5  The  bilinear  form  is  hounded  on  q  x  0,  be., 

Be{u,v)  <  C'||u||i,g||u||i_g-i,  for  all  u  e  HI  q,  v  G  Lbg-yo- 
Proof.  Let  u  G  q  and  v  G  q-  Then 


Bq{u,v)  =  /  [u'u' +  (1  +  0)uu]  dx 

Jr 

=  f  [g^^^u'g-^^^v'  +  (1  +  +  ey/^g-^/^]  dx 

Jr 

<  C[  f  {gu'"^  +  {I +  e)guydxY^'^[  f  {g~^v'^  +  {I  +  e)g~^vy  dx] 
Jm.  Jm. 


=  C||u||l,g||u||l,g-1.  □ 


Lemma  5.6  Suppose  ^  <  1.  Then  there  is  a  constant  C  >  0,  which  depends 
2 

on  such  that 


inf  sup 


Be{u,v) 


9  •*-  >© 


>  C  >  0. 
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Proof.  Suppose  u  €  Hg  q.  We  consider  v  =  gu.  Now, 
Bq{u,v)  =  (  [u' v'  +  Quv]  dx 


[u' {gu'  +  g'u)  +  Qgu^]  dx 
[gu'"^  +  Qgu^]  dx  +  f  uu' g'  dx. 


(5.35) 


Now,  for  e  >  0, 


/  uu' g' dx\  =  I  /  guu' {  —  )  dx\ 

Jr  Jr  9 

^  <^  [  \g^^‘^ug^^'^u'\dx 


<  a 

and,  therefore  from  (5.35),  we  get 

.'2 


[e  f  gu'^  dx -\ —  f  gu^  dx], 

Jr  c  Jr 


Bq{u,v)  >  I  {gu'^  +  Qgu^)  dx  —  a[e  I  gu'^  dx -\ —  /  gu^  dx] 

Jr  Jr  £  Jr 

=  (1  —  ae)  f  gu'^  dx  +  {1 - p^)  f  Qgu^  dx.  (5.36) 

Jr  Jr 

We  choose  e  such  that  ae  <  1  and  <  1,  and  therefore  from  (5.36),  we  have 


Be{u,v)  >  Cillulli 


where 


9’ 

a 


Cl  =  min[(l  -  ae),  (1 - p-)]  >  0. 

eB 


(5.37) 

(5.38) 


We  next  show  that  ||u||i^g-i  <  C'2||'u||i^g.  First  note  that 
j  g~^v'^  dx  =  j  g~^{gu'  +  g'uy  dx 

Jr  Jr 

=  I  gu'^  dx  +  I  g~^g'^u^  dx  +  2  f  g'uu'dx.  (5.39) 
Jr  Jr  Jr 


Now, 


and 


f  g  ^g'^u^dx=  f  g{—)^u^dx<a^  (  gu^  dx,  (5.40) 

J]R  J'R  9  Jm. 

2  /  g'uu'dx  =  2  f  g{  —  )uu'dx  <  2  f  ]ag^^'^ug^^'^u']dx 

Jr  Jr  9  Jr 

<  I  {gu^  +  a^gu^)  dx.  (5.41) 
Jr 
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Therefore  using  (5.40)  and  (5.41)  in  (5.39),  we  get 


g  dx  <  2  f  gu'^  dx  H — ^  j  Qgu^  dx. 
;  Jv:  c)  /r 


(5.42) 


Thus,  combining 

0  [  g~^v^  dx  =  Q  [  g~^g^v?dx  =  Q  [  gu^  dx 


with  (5.42),  we  get 


||v||i„-i  <  2  /  gu''^  dx  +  {1  + ‘^^)  [  Qgv?  dx 

./ra  ^  ./TO 


Since  ^  <  1,  therefore,  from  (5.43)  we  have 

Thus,  V  e  e>  combining  (5.37)  and  (5.44)  we  get 


.  .  Be{u,v) 

inf  sup  j— - 

9  ^  ,Q 


>C>0 


where 


min[(l-ae),(l-  ^)]  ^ 

73 


(5.43) 


(5.44) 


We  now  prove  the  inf-sup  condition  on  Sh  x  Sh-  In  the  proof,  we  will  use 
the  function  di{x),x  G  /((  and  i  G  to  denote  the  following  similar  functions: 

gi  -  g{x)  9i+^  -  9{x)  gi+i  - 


\/9H9i^  ^54+1 5(7 

where  Ik  G  A^.  It  is  easily  seen  from  the  definition  of  g{x)  that 

\di{x)\  <  Cah 


(5.45) 


2 

Lemma  5.7  Suppose  ^  <  Ci  and  ah  <  C2,  where  Ci,C2  are  sufficiently 

small.  Then  there  is  a  constant  C  >  0,  independent  of  u,  v,  and  h,  hut  may 
2 

depend  on  k  and  such  that  for  h  small  enough, 


inf  sup 

«eSh  y^Sh 


Be{u,v) 

l|u||l.gl7lll,S-i 


>  C  >  0. 


(5.46) 
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Proof.  Let  u  =  ll^^lli.g  ^  oo-  Then  for  x  S  I^,  we 

have  u  =  Z^jeAj  Since  Z^ieA^  4>i' (x)  =  0  for  a;  G  /^,  we  have 

«'(a^)  =  H  (x)  =  '^{df  -  cfj^f  (x),  a;  G 


ieA^ 


iGA^ 


where  Ik  G  is  a  fixed  integer  for  given  k. 

We  now  choose  v  =  J^iez'^idi+^di  S'/I,  and  as  before,  for  x  G 

*GA^ 

=  -  c^,gi,+  i)(p-'(x) 

^eA^ 

=  J2  +  c’f^  (ft+i  -  gi,+L)di\x). 

iGAJ 


»ga^ 


Now, 


j  uv'dx=  j  gu^  dx  +  j  u{v'  —  gu)dx.  (5-47) 

JlR  t/M  JlR 


For  e  >  0,  we  have 

/  t6'(w'  —  (/m')  dx 

Jm. 


<  el  gu'^  dx  -\ — 

./to  e 


1  /■  (t''  -  gu'f 


■  dx.  (5.48) 


Now,  from  the  definition  of  v'  and  u' , 

1 


-{v'  —  gu'Y  dx  = 

hi  g 


jh 

-'fc  iGA' 


E(4-4) 


5i+i  g 


g 


1/2 


<  S  "-A 

ieAt  ^ 


<  C  [J2{c>l-ct) 


3^+h  9  ih'-\2 


I  jh 

T  ieAi; 


jl/2 


(jii  Y  dx 


E  ^A±-^^YYdx.  (5.49) 


jGAJ 


The  first  term  of  the  RHS  of  the  above  inequality,  employing  (5.45)  and  (5.20), 
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/ 

Ic  i£A’i 


/  T^ 


-  ^  f  ^tf9^+^dj{(j)ff  dx 

<  Ca^/i^  ^  (cf  -  /  (</’f  )^ 

<  Ca2/i2l  ^ 

idAl 

<  Ca^h^l  Y.  (4i-c'‘)'ft+i, 

i.(i+i)eA^ 

where  C  is  independent  of  a,  h,  but  depends  on  k. 

The  second  term  of  the  RHS  of  (5.49),  employing  (5.45),  gives 

/  IS 


_  t  h\2  /  ^*+5  ^'fc  +  5  //i'l2  , 

'*  ■+*4£(94  +  5)'''V'"‘^-' 

<  Y  [  dj{(l)^'fdx 

Ah  J  lu 


i&Ai  k. 


<  Ca^h{clfgi,,  (5.51) 

where  C  depends  on  k,  but  is  independent  of  a,  h.  Therefore,  from  (5.49), 
(5.50),  and  (5.51)  we  have 

[-{v'-gu')^dx  <  Y  i.cY  -  c'lf  9h+\ 

d  i,p+i)gAj 

+Ca^h{c'^J^gi^. 

Now  summing  the  above  inequality  over  k  €  1a,  and  using  (5.27)  and  (5.28),  we 
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get 


[  -{v' —  gu')"^  dx  =  [  -{v' —  gu')^  dx 

^  fcez  ^ 

fcez  i_(i+i)gA^ 
feez  igAj 


< 


iez 


+C'a^  E  /,  dx 


<  Ca^h?  I  gu'^  dx  +  Ca‘‘  I  gu"^  dx. 


kez''^k 
„/2 


(5.52) 


Then  from  (5.48)  and  (5.52)  we  have 

/  u'{v'  —  gu')dx  <  e  gu'^  dx 

Jr 


+  -[Ca^h^  f  gu'^dx  +  Ca^  f  gu^  dx] 

e  Jk  Jm 


=  ie  + 


Ca^h 


2u2 


)  /  gu'^  dx  + 


Co? 


Qgv?  dx.  (5.53) 


We  next  consider 

0  f  uvdx  =  Q  (  gv?  dx  +  Q  f  u{v  —  gu)  dx. 


(5.54) 


For  ei  >  0,  we  have 

u{v  —  gu)  dx 


1  9  '  , 

<  ei  [  gu'^dx+-  [  ~  dx.  (5.55) 

Jk  ^1  Jm  9 
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Now, 


r  {v-  guf 

I  It  9 


dx  = 


^  ieAj 

-  ^  [  Y1  (<^'^y9rd'i4>f  dx 

^  I  u  .  _  ^  h 


fc 


<  Ca^h'^h  ^  (c'*)^5i. 


Therefore  using  (5.27),  we  get 
f  {v  -  guf 

Jr  9 


dx  ^  YJ 


{v  -  guf 


dx 


kez-'^k  d 


k^lj 


idAl 


<  Ca^h?  I  gu‘‘  dx. 


Thus  from  (5.55),  (5.56),  we  have 


0|  f  u{v  —  gu)  dx\  <  {ei -\ - )  f  Qgv?  dx, 

JR  ei  JR 

and  combining  (5.47),  (5.53),  (5.54),  and  (5.57),  we  get 
\Bq{u,v)\  >  [  gu^dx  +  Q  [  gv?  dx 


—  I  /  u  {v  —  gu)  dx\  —  Q\  /  u{v  —  gu)dx\ 


>  (1-e- 


Ca'^h’^ 


)  j  gu'"^  dx 

JR 


M  Ca2/i2  Ca\  f  ^  2.. 

+(l-ei - ^)  /  Qgu^dx. 

El  £0  Jr 

Now  we  can  choose  e  and  ci,  for  sufficiently  small  h,  such  that 

\Be{u,v)\  >  Cillwll^^g, 

2 

where  Ci  >  0,  since  ^  <<  1,  ah  «  1  by  assumption. 


(5.56) 


(5.57) 


(5.58) 
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We  now  show  that  ||u||i,g-i  <  CHuHi^g.  From  the  definition  of  v' ,  we  have 

i&A'i 

+cik  gi^  +  l)(t)’i'f  dx 


i&A'i 


< 


'= 

(5.59) 


k  i^Ai 


Now, 


k  *eAj  ^ 


E  (cf  -  <)9''V?'  +  z  (e?  -  ^  0?? 

<  C*  /  +  ~  ^h)^9^+^dU(t>i')^  dx 

^k  iGAk  ieAk 

<  C  gu'^  dx  +  C  {c^  -  c^J'^g^+^dj{4>^'f  dx.  (5.60) 


k  i^Ak 

Also  using  (5.45)  and  (5.20),  we  have 


Y  i  /  d*  ((/''‘V  dx  <  Ca^h'^^  Y  “  'AJ^9^+^ 

leAk  ^k  i^Ak 


<Ca^h^^  Y  (cti-cf)Vi-(5-61) 

i,i+l&Ak 

Therefore,  using  (5.60),  (5.61)  and  (5.51)  in  (5.59),  we  get 


hi 


g  ^v'"^  dx 


-  ^  f  ^9u'^  dx  +  Ca‘^h{cfj‘^gi^+Ca^h‘^^  Y  {cY  -  <a)'^9z+^ 

^k  i,j+ieAj 

<  C  1^9'^'  dx  +  Ca'^h  "Y^  {c^)^gi  +  Ca^h^  —  'Y  (cf+i  —  Ci*)^5i+i 

A  i,i+leAk 

<  C  [  gu'^dx  +  Ca^  f  gu^dx  +  Ca^h^^  "Y  (c(Yi  ~  c(‘)^5j_i_  i . 

Jit  Jit  d  2 


i.i+leAj 
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Now,  summing  the  above  inequality  for  all  k  and  using  (5.20),  we  get 

j  g~^v'^  dx  <  C{\  +  a^h^)  j  gu^dx  +  Ca^  j  gv?  dx.  (5.62) 

jIR 


Again, 


g  ^v'^dx  =  f  g  ^C^c^gicp^f  dx  =  f  c’l dx .  (5.63) 


Now  using  (5.45),  we  get 

,1/2 


^  i^Ai 


I  Ip- 


h9i  ~  9-/112 


leAt 


9 


1/2 


(j)^]  dx 


-  E  +  C  /  ^  (cf)^9* 

<  C  [  gv?  dx  +  C  {c^Y gi  [  d‘f(j)'l‘^  dx 

<  C  f  gu^dx  +  Ca^h^h  ^  {c'lYgi 

^k  l^A\ 

<  C  [  gv?  dx  +  Ca^h^C  [  gu^  dx 

Jii  Jii 

<  C{l  +  a^hY  [  gu^  dx, 

Jit 

and  therefore,  from  (5.63)  and  the  above  inequality. 


(j)^  dx 


^V^dx  <  C'E  /  (E  '=^^72^*'“)"'^^ 


/  Th 


<  C{l  +  a^hY'^2i  [  911^  dx 
kezJit 

=  C{l  +  a^hY  [  gu^  dx. 

Jr 
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Thus  combining  (5.62)  and  above  inequality,  we  have 

lkll?,g-i  =  [  dx  +  e  f  g~^v^dx 

JR  JR 

<  C{1  +  a^h^)  (  gu'^dx  +  Ca^  f  gv?  dx 

JR  JR 

a^h^)  f  Ogu^  dx 

JR 

<  C{1  +  a^h^)  f  gu'^  dx 

JR 

2  p 

+  — \- C{1  +  a^h^)]  /  Qgu^  dx 

“  JR 

<  c(l  +  +  ^)\\u\\l^  <  CMlr  (5-64) 

Finally,  combining  (5.58)  and  (5.64)  we  get  the  desired  result.  □ 

Projection  with  respect  to  Bq{u,v): 

Suppose  u  G  Hg  Q  and  let  Pqu  be  the  projection  of  u  onto  Sh  defined  by 

Bq{Pqu,v)  =  Bq{u,v),  for  all  v  G  Sh- 

The  projection  Pqu  exists  (see  [11]),  and  it  is  clear  from  Lemmas  5.7  and  5.5 
that 

||-PeM||i.3  <  c  sup  <  C'||M||i,g.  (5.65) 

We  first  note  that  for  fixed  h,  a,  and  0,  the  polynomials  belong  to  the  space 
^g,e-  Moreover,  for  fixed  ft.,  a,  and  0,  we  can  also  show,  using  (5.27)  and 
Remark  5.6,  that  ih(x^~^^)  G  q,  where  X?i(a;*+^)  is  the  interpolant  of  x^^^, 
as  defined  in  (4.35). 

We  now  present  some  simple  facts  about  polynomials  and  periodic  functions. 
Lemma  5.8  Let  the  shape  functions  {(p'fjiizz  be  reproducing  of  order  k.  Then 

(a)  Pqx’'  =  x*,  0  <i  <  k  (5.66) 

(b)  Peih{x'^+^)=ih{x'^+^),  (5.67) 

where  Ih{x^^^)  is  the  interpolant  of  x^^^  as  defined  in  Section  j. 

The  proofs  of  these  facts  are  immediate. 

Lemma  5.9  Suppose  f  €  q  is  periodic,  i.e.,  f{x  +  x\)  =  f{x)  for  all  k. 
Then  P^f  is  also  periodic. 

Proof  Let  f{x)  =  f{x  +  x^).  Then  [PefKx)  =  [PefKx  +  x^).  Now  f{x)  = 
f{x)  since  /  is  periodic,  and  thus  from  the  uniqueness  of  the  the  projection  Pq, 
we  have  [PQf]{x  +  x^)  =  [Pe/](a;),  i.e.,  Pef  is  periodic.  □ 
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Remark  5.8  We  note  that  if 

iez 

is  a  periodic  function,  i.e.,  v{x  +  x^)  =  v{x)  for  any  k,  then  v  is  a  constant. 
This  could  be  shown  as  follows:  Since  v{x  +  x\)  =  v{x),  we  have 

v{x  +  ^  ^  ^  cf(/)^(a;)  =  v{x), 

iez  iez  iez 

which  implies  that 

^[c'Vfc  -  cf]^f(a;)  =  0,  for  all  a;  gM. 
iez 

Using  (5.19),  we  can  show  that  {4>i}i^i.  are  linearly  independent  in  M.  Thus  we 
infer  from  above  that  =  C  (constant),  for  all  i  G  Z.  Recalling  that 

form  a  partition  of  unity,  we  get  v{x)  =  =  C. 

We  now  define 

ef+i(x)^x'=+i-Pex'=+^  (5.68) 

^®+i(a^)  and  the  next  result  will  play  a  central  role  in  the  final  result  of  this 
section. 

Lemma  5.10  Let  ^f_^i{x)  be  as  defined  in  (5.68)  and  consider  (fj:_^i{x)  =  x^^^  — 
defined  in  (4.36).  Then 

^®+i'(a:)  =  Cfc+/(a:).  (5.69) 

Proof.  We  first  note,  from  the  definition  of  ,  that  £/^j^i{x)  =  x^^^  — 

XhX^^^.  Now,  using  (5.67),  we  have 

it+i  =  -  Pex'^+^ 

=  x’^+^  -  +  ihx’^+^  -  Pex^+^ 

=  ek+i-Pe[x'^^^-ihx'^^^] 

=  4\i-^eKti]-  (5.70) 

But  we  know  from  Lemma  4.1  that  £)(._^_i{x)  is  periodic,  and  therefore  from 
Lemma  5.9  and  Remark  5.8  we  infer  that  Pei'Cfc+i]  is  a  constant.  Thus,  from 
(5.70),  we  get 

?®+i'(a:)  =Cfc+/(a:), 

which  is  the  desired  result.  □ 

Proof  of  Theorem  5.1: 

The  proof  will  be  given  in  several  steps. 
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1.  Let  E  be  the  extension  operator  satisfying  (3.51).  Then,  for  x  G  Bh  = 
Bh{0),  we  have 

[uo-Uh]{x)  =  [uo-Pe{Euo)-{Euh-Pe{Euh)}  +  {Pe{Euo)-Pe{Euh)}]{x), 


and  therefore, 

(uq  -  Uh'){x)  =  {mo  -  Pe{Euo)y{x)  -  5'h{x)  +  p'h{x),  (5.71) 

where 

Sh  =  Euh-Pe{Euh);  (5.72) 

Ph  =  Pe{Euo)  -  Pe{Euh).  (5.73) 

Since  uq  =  Euq  in  Bh{0),  from  Taylor’s  Theorem  we  have 

Euo{x)  =  .  x^  +  +  Rk+i{Euo){x),  (5.74) 

^  j!  (k  +  iy. 

where  Rk+i{Euo){x)  is  the  remainder  given  by 

Rk+i{Euo)ix)  =  Hx  -  t)>^+\Euoy'^+^\t)  dt.  (5.75) 

(fc  +  1)!  Jo 

Since  Pq  is  a  linear  operator,  we  have 


Pe{Euo){x)  =  Y  "^^Pex^  +  7,  ^  J,  Pex’^^^  +  PeRk+i (Euo) (x) .  (5.76) 
^  j!  (k+l)l 

We  know  from  (5.66)  that  Pqx^  =  x^O  <  j  <  k.  Therefore,  by  first  subtracting 
(5.76)  from  (5.74),  then  differentiating  the  identity,  and  finally  using  (5.69),  we 
have 

{Euq  -  Pe{Euo)y{x) 

('Q'l 

=  -  Pex'^+^y{x)  +  [Rk+i{Euo)nx)  -  [PeRk+i{Euo)nx) 

(fc+1)  /Q\ 

=  ""(1+1)!  +  [Rk+i{Euo)nx)  -  [PeRk+i{Euo)nx) 

=  (fc+1),  (^)  +  [Rk+i{Euo)nx)  -  [PeRk+i{Euo)nx).  (5.77) 

Thus  from  (5.71),  (5.77),  and  using  eh{x)  =  [uq  —  Uh]{x),  we  get  for  x  G  Bh, 

(fc  +  1)  /^\ 

=  [Rk+i{Euo)]'{x)-[PeRk+i{Euo)nx)-6'k  +  pl  (5.78) 
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2.  From  (5.75),  we  have 


[Rk+iiEuo)nx)  =  i  J\x  -  dt, 

and  since,  I|wo||^'=+2(52h)  —  have  for  x  €  B2H, 

f  \[Rk+i{Euo)]'\'^  dx  <  [  g\[Rk+ii.EuQ)]'\'^  dx 

J  Bh  Boh 


'  B2H 


(5.79) 


'  dx 


Similarly,  again  from  (5.75),  we  get 

[  \Rk+i{Euo)f  dx  <  [  g\Rk+i{Euo)\ 

J  Bh  B2H 

3.  It  can  be  shown  from  the  definition  of  g{x)  that,  for  0  <  j  <  /c  +  1, 


(5.80) 


gx 


^^dx=  /  dx  < 


(5.81) 


I2H 


f2H 


where  C  depends  on  A:  +  1.  Now,  from  (5.65)  we  get 

\[PeRk+i{Euo)]^^  dx  <\\PeRk+i{Eu^)\\lg  <C\\Rk+i{Euo)\\lg.  (5.82) 


I  Bh 


We  note  that,  from  (5.74),  we  have 


[Rk+i{Euo)nx)  =  {EuoYix)  -  ^ 


i=o 


(j  +  1)! 


Therefore,  using  (5.81)  and  the  fact  that 

I  ^OO 

gKEuoYY  dx  <  e~°'^  /  \{EuoYY  dx  <  e~°‘^\Euo\Hi(K), 

J2H 


f2H 


we  have 


f2H 


g\[Rk-\.i{Euo)y\  dx 


(5.83) 
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where  C  depends  on  k.  Similarly,  we  can  show  that 

/  9\[Rk+i{EuQ)]'\^  dx  <  Ce  “^{||Mollffi(n)  +  C'||'Uo||^M-2^g^^^}, 

J  — OO 

which,  together  with  (5.83)  imply  that 
f  g\[Rk+iiEuo)]'\^dx  <  Ce-“^{||wo|ll^i(o)  +  (5.84) 

Using  similar  arguments,  we  can  show  that 

[  g\Rk+i{Euo)\^dx<Ce-°‘^{\\uo\\l^^^)  +  C\\uof,,+2,  A.  (5.85) 

JR-B2h  00  V 

Now  combining  (5.79),  (5.80),  (5.82),  (5.84),  and  (5.85)  we  get 


/  \[PeRk+i{Euo)y\^dx 

J  Bh 

<  C\\Rk+i{Euo)\\lg 

=  C  [  gllRk+iiEuo)]^"^  dx  +  Ce  /  5|i?fe+i(UMo)P  dx 
Jr  Jr 

+C'(l  +  0)e  “^{||uo||^i(0)  +  C'||Mo||^^+2^g^^j}.  (5.86) 

4.  We  first  note  from  (5.72),  that 

f  Sy^dx<\\6alg  =  \\Euf,-Pe{Euh)\\lg.  (5.87) 

JBh 

Let  PQ{Euh)  =  i^Eufi  +  £■  Then  S  €  Sh-  Now  from  Lemma  5.5  and  the 
definition  of  Pe,  we  have  for  all  v  G  Sh, 

Bq{£,v)  =  Be{Ps{Euh)  -ihEuh,v) 

=  Be{Euh -ilEuh,v) 

<  C\\Euh-i*hEuh\\iJv\W^g-., 

and  hence  from  Lemma  5.6,  we  get 

Bq{E,v) 


Thus, 


lli^lli.g  <  C  sup  <  C\\Euh-ThEuh\\i,g. 


\\EUh  -  Pe{EUh)\\l,g  <  \\EUh-i*hEUh\\l,g+\\£\\l,g 

<  C\\Euh-ThEuh\\i,g.  (5.88) 
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We  now  estimate  the  RHS  of  the  above  inequality.  We  first  note  that  Euh{x)  = 
Uh{x)  for  X  G  ^1.  Consider  £2  C  such  that  (see  Remark  5.4) 

B2h  C  0  and  i^Euh\n  =  Euh\n  =  Uh\n. 

Therefore,  from  Lemma  5.2  and  using  (5.5), 


g[{Euh  -  ilEuh)'f  dx 


g[{Euh  -  ilEuh)']'^  dx 


Jr  Jr—q 

<  Ce  °‘^[\Euh\jji(j^-)  +  j^\\Euh\\‘i^(j^-)] 

-  pe-“'^||uo||^i(n).  (5.89) 

Similarly,  we  can  show  using  Remark  5.5  that 

[  g[Euh-ilEuhf  dx  <  C'e"“'^||uo||i2(n), 

Jr 


and  thus  combining  it  with  (5.89),  we  get 

CQ 

\\Euh-l*hEuh\\l^g  <  -^e'  —  lluolli^qo). 
Now  from  (5.87),  (5.88),  and  above,  we  get 

C© 


IBh 


S'Jdx<'^e--^ 


2 

H^Q)- 


5.  We  first  note  from  (5.73)  that 

[  pjdx  <  WphWlg  =  \\Pe{Euo)  -  Pe{Euh)\\l, 


IBh{0) 

Now  using  (5.4),  we  have  for  all  v  G  Sh, 


(5.90) 


(5.91) 


Be{ph,v) 

=  Bq{PqEuq  -  PeEuhjv) 

=  Bq{Euo  -  Euh,v) 

=  B^{Euo  —  Euh,  v)  +  B^~^ {Euo  —  Euh,  v)  +  QD^{Euo  —  Euh,  v) 

=  B^{uo  —  Uh,v)  +  B^~^{Euo  —  Euh,v)  +  QD^{Euo  —  Euh,v) 

=  B^-^{Euo  -  Euh,v)  +eD^-^{Euo  -  Euh,v)  +eD^{Euo  -  Euh,v) 
=  Bl-^{Euo-Euh,v)+eD^{uo-Uh,v).  (5.92) 
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Also,  for  V  G  Sfi 


Bq  ^{Euo  -  Euh,v) 


i-n 


[{Euq  —  EuhYv'  +  Q{Euo  —  Euh)v]  dx 


<  C\\\Euo  —  |||u|||i,g-i,K-n 

<  C\\\Euo  -  Iklll,g-1, 


(5.93) 


where 


ll,g-i.R-f2 


g  dx  +  Q 


g  dx] 


\\\Euo  -  Euh\\\Yg^R_Q  =  (  giEuo-Euh)' 

JR—Q 


l-Q 


dx 


+0  /  g{Euo  —  Euh)^  dx. 


i-a 


From  the  definition  of  g{x),  we  can  show  that 

\\\Euo  —  Euh\\\Yg^R-ii  <  e  °‘^Q\\Euo  —  Euh\\‘Hi^B^_Q) 

<  Ce-^^muo-UhlYrn^n) 

<  Ce-“^0||uo||?,.(a). 

Now  using  the  definition  of  g{x)  and  (5.6)  with  R  =  2H,  we  get 


(5.94) 


/  g{uo-Uhfdx  =  /  g(uQ-UhYdx+  /  g(uQ  -  UhY  dx 

J  Vl  B2H 

—  ocH  I 


and  therefore, 


0 


0 


<  ||mo  —  +  e  “  ||mo  —  M?i|lL2(a) 

-D^{uo  -  Uh,v) 


< 


l^'lll.g-i  Jn 

0 


(uo  —  Uh)v  dx 


Plli,g-1  \Jn 


\  1/2 


g{uo-Uh)  dxj  IJ  g  dx 


1/2 


<  e-^Ch'^+^H-Y\uo\\H^+i^a)  +  e-^e-^^/^uo\\HHn)-  (5-95) 

Now,  from  the  inf-sup  condition  (5.46)  and  using  (5.92)  and  (5.93),  we  have 

Be{ph,v) 


\\Ph\\l,g  <  CSUP 

v&Sh  ll^'lll.g-i 

<  C\\\Euq  -  Euh\\\i,g,R-n  +  sup 


0 


eSh  ll^lll.g-i 


-D^{uo  -  Uh,v), 
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and  thus,  using  (5.91),  (5.94)  and  (5.95),  we  have 

/  P'ndx  <  WpuWI, 

JBh 

<  (5.96) 

6.  We  first  note  from  (5.69)  that  {x)  where  is  defined 

(fc+1)  /p.-. 

in  (4.36).  Let  T{uo)  =  -p^ipYyp.  Then  from  (5.78),  we  have 
eh'(a;)  -  T(uo)Cfc+/(a;) 

=  [Rk+i{Euo)]' (x)  -  [PeRk+i{Euo)]'{x)  -S'f^  +  p^, 
and  therefore,  from  (5.79),  (5.86),  (5.90),  and  (5.96),  we  have 

[  (eh'-T{uo)^^+A\x 

JB„  ^  ^ 

<  c[  \[Rk+iiEuo)]'\^dx  +  C  [  \[PeRk+iiEuo)]'\‘^dx 

J Bm  'J  Bm 


>Bh 

+C  f  6'u^dx  +  C 

Jb„ 


Phdx 


r'A 

+C'(l  +  0)e-“^{||uolll,i(n)  +  ^ 

+Ce-“^0||uo||?^i(o)  +  C'0h2'=+2i7||wo||^.+,(f,) 

<  C  +  (1  +  0i72)H2fc+2^  (J^  0)e-“'^ 

0 


2 

//1(a) 


+  pe-“^  +  M{uo),  (5.97) 


where 


M{uo)  —  ||■uo||^fc+l(Q)  + 

We  will  now  choose  a,  0,  and  R,  where  E[  =  and  7  <  1.  First  we  choose 
7  such  that 

jjk  +  2  ^  J^k  +  l^  (-5  gg^ 

which  implies  that 

/j7(fc+2)  ^  J^k+l 

or,  7(fc  +  2)  =  fc+1 

fc+1 
fc  +  2 


or,  7  = 


<  1. 


Let  e  >  0,  which  depends  on  7,  be  such  that  e*  =  l  —  7  —  e>0.  We  will  now 
choose  a  such  that 


e““^  <  =  /i2fe+4+37+2e 


(5.99) 
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This  implies  that 


a  >  Ci{\Tih-^)h-'^ , 

where  Ci  =  2fc  +  4  +  87  +  2e.  Since  h~'^  >  In  for  small  h,  we  take 

a  =  Cih-^^+‘^\  (5.100) 

We  now  choose 

0  =  C2>Ci.  (5.101) 

We  note  from  (5.100)  that  ah  =  =  Cih^  <  1  for  small  h,  and 

lim/j^o  cth  =  0.  Thus  ah  can  be  made  sufficiently  small;  this  was  one  of  the 
assumptions  in  Lemma  5.7.  Also  by  choosing  C2  large  enough  in  (5.101),  we 
can  make  ~  =  (C'i/C'2)^  <<  1,  ie.,  sufficiently  small,  which  was  another 
assumption  in  Lemma  5.7.  Thus  the  conclusion  of  Lemma  5.7  is  true  for  the 
choices  of  a  and  0  given  in  (5.100)  and  (5.101),  respectively. 

Now,  for  these  choices  of  7,  a,  and  0,  we  have 


Qh‘^k+2  ^ 

0/12(7+0/12/0^2(1-7-6)  ^  fj2^2k+2e*  ^ 

(5.102) 

Using  (5.98)  and  (5.102)  we  have 

0^2^2fe+2 

=  0i72fc+4  =  0/i2fe+2  ^  pj2j^2k+2e-  ^ 

(5.103) 

Also  from  (5.99),  we  get 

0e-“-f^ 

<  /l2'=+4i70/i27+26  <  (j2f^2k+2jj^ 

(5.104) 

and 

0  -aH 

<  /l2'=+2i70/i27+26  =  (j2j^2k+2^_ 

(5.105) 

Thus,  using  (5.98),(5.102)-(5.105)  in  (5.97),  we  get 

W  -  T(uo)5ti'llL.(B«)  <  Ch'^+^'H^/^M{uo)K  (5.106) 

and  hence  using  (5.7)  with  p  =  H,  we  have 

||e/i'  -  r(Mo)^fc+l  \\L2iBH)  ^  p.e* 

l|e.'||L.(B„) 


where  M('Uo) ^ /||uo||//fc+i(o)  <  C,  which  is  the  desired  result.  □ 

Remark  5.9  The  balancing  of  various  terms  in  item  no.  6  of  the  proof  of  The¬ 
orem  5.1  is  similar  to  balancing  used  in  the  superconvergence  of  FEM  solution 
(see  [19]. 
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Remark  5.10  Assuming  that  our  superconvergence  result  is  valid  in  Loa,  i.e., 
assuming  that  for  x  G  Bm  there  exists  e*  >  0,  such  that 


we  see  that  the  zeros  of  ^fc+i^(f)  are  the  superconvergence  points.  In  Figure 
5.1, 


Figure  5.1:  The  plot  of  ^2  (2/),  0  <  y  <  1  for  (a)  RKP  shape  functions, 
reproducing  of  order  fc  =  1,  corresponding  to  the  conical  weight  function 
with  1  =  2,  R  =  l.S  (b)  standard  “tent”  functions  used  in  FEM. 

we  have  presented  the  plot  of  ^k+i'iv)  the  RKP  shape  functions,  reproducing 
of  order  fc  =  1,  with  respect  to  the  weight  function  w{x)  given  by  (4.4)  with  I  =  2 
in  1-d.  We  have  also  included  the  plot  of  ^k+i'{y)j  k  =  1  (the  dashed  curve)  for 
the  standard  tent  functions  that  are  used  as  shape  in  FEM.  We  note  that  ^2'{y) 
for  the  tent  function  has  only  one  zero,  where  as  ^2  {y)  has  5  zeros.  Thus  the 
superconvergence  points,  for  the  RKP  shape  function  could  be  distributed  quite 
differently  than  the  corresponding  points  for  standard  tent  functions  in  FEM. 


6  The  Generalized  Finite  Element  Method 

The  idea  of  the  Generalized  Finite  Element  Method  (GFEM)  was  first  intro¬ 
duced  in  [16]  to  address  elliptic  problems  with  rough  coefficients.  This  idea  was 
later  extended,  and  called  the  Partition  of  Unity  Method  (PUM)  in  [17]  and 
[62].  In  the  current  literature,  PUM  is  referred  to  as  GFEM  ([82], [83]).  In  this 
section,  we  will  first  describe  the  GFEM  and  present  the  relevant  approxima¬ 
tion  results.  We  will  then  discuss  the  selection  of  an  optimal  or  near  optimal 
approximating  space,  to  be  used  in  the  GFEM,  in  certain  situations. 
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6.1  Description  of  GFEM  and  Related  Approximation 
Results 


In  this  section  we  will  discuss  the  GFEM  in  the  context  of  general  particle-shape 
function  systems,  which  were  discussed  in  Section  3.3.  Suppose  uq  is  the  solution 
of  our  model  problem  (2.1),  (2.2)  (or  (2.3)).  We  consider  a  family  of 

particle  shape  function  systems  satisfying  assumptions  A1-A7  with  fc  =  0  and 
=  /;  assumption  A5  then  reads 


4>x{^)  =  1;  for  all  X  G  M".  (6.1) 

x&X" 


The  partition  of  unity  (6.1)  is  the  starting  point  of  GFEM.  We  will  need  addi¬ 
tional  assumptions  on  {Af'^}i/gAr,  namely. 


Wx 


.)  <  Cl 


(6.2) 


and 


< 


C2 

diam(r7^)  ’ 


(6.3) 


for  all  X  G  X'' ,  and  all  v  G  N .  In  (6.3),  we  implicitly  assume  that  q  >  n/2.  We 
also  assume  that  there  is  a  constant  C  such  that 


dia,m{rix)  <  C,  for  all  x  G  X''  and  for  all  v. 

For  each  x  G  X'',  we  assume  that  we  have  a  finite  dimensional  space 
of  functions  that  has  good  approximation  properties.  We  refer  to  as  local 
approximating  spaces.  We  define  a  set  of  particles  Aq  ,  namely, 

A^  =  {xGX^:y^nn^0},  (6.4) 

for  each  v  G  N .  From  (6.1)  we  have 

=  1,  for  all  X  G  fl.  (6.5) 

For  an  approximating  space  on  we  then  consider 

=  {^In  ■  ^  =  51  ^  (6-6) 

The  GFEM  is  the  Galerkin  method  (2.7)  with  B  =  B  and  S  =  and 
we  will  denote  the  approximate  solution  us,  obtained  from  GFEM,  by  ugfem- 
When  GFEM  is  used  to  approximate  the  solution  uq  of  the  Neumann  problem, 
Vx  can  be  any  finite  dimensional  subspace  of  But,  when  GFEM  is  used 

to  approximate  the  solution  uq  of  the  Dirichlet  problem,  with  the  boundary 
condition  (2.3),  the  functions  in  Vx  are  required  to  satisfy  =  0,  for 

particles  x  for  which  \f]x  H  n|  >0.  Thus  the  approximating  space  V'  C  i7o(n). 

Our  next  theorem  states  an  approximation  result  for  V'^ .  We  will  follow  the 
ideas  presented  in  [15,  16,  17,  62,  82,  83]. 
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Theorem  6.1  Suppose  u  G  and  suppose,  for  all  x  G  Aq,  there  exists 

ipf  G  Vf!  such  that 


Then  the  function 


satisfies 


||m 


\\u- 'ipxWL^iriinn)  <  ei(x), 
\\^{u-ifx)\\L2{vinQ)  <  <^2{x). 

«ap=  E 

Uap\\L2m  ^  ' 


(6.7) 

(6.8) 

(6.9) 


(6.10) 


and 

||V(w 


U, 


ap)\\L2{^) 


1/2 


(6.11) 


Proof.  We  will  prove  only  (6.11),  since  (6.10)  can  be  proved  similarly. 
Since  (ff,  for  x  G  form  a  partition  of  unity  for  fl  (see  (6.5)),  we  have 


II  V(m  —  Uap)||i2(f2) 

=  ||V  ^  C(«-V’.)llL(o) 

<2||  (w-V'.)VC||i,(0)+2||  E  CV(zi-V'.)llL(n)-(6.12) 


For  any  x  G  fl,  the  sums  ~  4‘x^{u  —  f’x)  have  at 

most  K  non-zero  terms  (see  Remark  3.4  and  (3.61)).  Therefore, 

I  ^  (u-  V'2;)V(/)^P  <  K 

x&A^  xeA'^ 

and 

I  Y  CV(u-^x)P<«;  Y 

xeA^  xeA^ 

Hence,  from  (6.12),  (6.7),  (6.8),  recalling  that  supp(^5()  =  pf.,  we  have 


|lV(w  —  Uap)\\L2(Q) 

<2kY  \\{u-i^x)^rjUn)  +  2^  Y 

x&A'^  x&A''^ 

=  2kY  ll(«-V’^)VCIlL(onp£)+2«  E  IICV(«-V'^)||L(anpi) 

x&A'^  x&A^ 


<2nY 


c. 


'diam(77^) 


^f4(x)  +  C^eUx)'), 
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which  is  the  desired  result.  □ 


Remark  6.1  We  note  that  ei(x),  e2(x)  in  (6.7),  (6.8)  depend  on  the  parameter 

V. 


We  will  now  show  that  both  the  terms  of  the  estimate  (6.11)  are  of  the  same 
order  with  additional  assumptions  on  .  These  additional  assumptions  depend 
on  the  boundary  conditions  of  the  approximated  function. 

Theorem  6.2  Suppose  uq  G  is  the  solution  of  the  Neumann  problem 

(2.1),  (2.2),  and  suppose  there  exists  iff.  G  Vf,  x  G  Aq,  such  that  (6.7)  and 
(6.8)  are  satisfied.  Moreover,  assume  that  for  x  G  Aq,  the  space  Vf  contains 
constant  functions  and  that 

inf  Ik  -  M\L2(v-nn)  <  C  (diam(r]f))  ||  for  all  v  G  H^{r]f  n  ft), 

Atil^  —  — 

(6.13) 

where  C  is  independent  of  x  G  X''  and  v.  Then  there  exists  iff  G  Vff  so  that 
the  corresponding  function. 


Uap  =  Y.  e 


satisfies 

\\uo-Uap\\Hi{a)  <c(^Y  )  (6-14) 

where  C  is  independent  of  uq  and  v. 

Proof.  Let  iff  G  Vff,  x  G  Af^,  satisfy  (6.7)  and  (6.8).  Define  iff=iff  +  rf, 
where  rji  G  M  satisfies 

Iko  -  Clk2fenn)  =  inf  |ko  -iff-  A||L,(,,j;nn)-  (6.15) 

—  A  tiK  — 

Since  VJf  contains  constant  functions,  it  is  clear  that  ifx  G  14-  Also,  from  (6.15), 
(6.13)  with  V  =  Uq  —  iff,  and  (6.8),  we  have 

Ik  -  '^x\\L2{v^nn)  <  C  diam(77^)  |jV(u  -  kDIU^feno) 

<  C  diam(?7^)  £2(0;).  (6.16) 

Let  Uap  =  '/'a-k®-  Recall  that  (ff,  x  G  Aq,  is  a  partition  of  unity  for 

n.  Then,  following  the  arguments  in  the  proof  of  Theorem  6.1  and  using  (3.61), 
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(6.2),  we  can  show  that 


Ik-  Wap|li,(o)  =  II  Ck-kkllL(O) 

<  kY^  IIC(w-k£)llL(n) 

=  ^  lkK(^^  —  kkllLaCnn??^) 

<  c  ^  ||(«-4)|k^(^n.i),  (6.17) 

and  using  (6.16)  in  this  inequality,  we  get 

Ik  -  UapWYn)  ^  C'  H  (diam(?7^))k^(x).  (6.18) 

£6^si 

Again,  following  the  arguments  in  the  proof  of  Theorem  (6.1),  and  using  (6.2), 
(6.3),  we  can  show  that 


||V(u  —  Map)||i2(n) 

<2k  ^  ll(«-4)V4llL(nn,s)+2«  E  IICV(u-4)||L(nnc) 


(diam(r7^))2 


k-CI 


2 

L2(nn77|^) 


+C  Ilv(u-4)||i^(^n,^).  (6.19) 


By  first  noting  that  \7{u  —  ip^)  =  \7{u  —  ip^),  and  then  using  (6.16)  and  (6.8)  in 
the  above  inequality,  we  get 

llV(u  -  fiap)llL(o)  <  C-  E  ^2(31).  (6.20) 

Combining  this  with  (6.18)  we  get  (6.14),  where  we  used  that  diam(?7(()  <  C  for 
all  X  G  X‘'  and  for  all  v.  □ 


Theorem  6.3  Suppose  uq  G  Hq{^1)  is  the  solution  of  the  Dirichlet  problem 
(2.1),  (2.3),  and  suppose  V)f ,  x  G  Aq,  satisfy  the  following  assumptions: 

(a)  For  all  x  G  such  that  pf.  n  dFt  =  0,  Vf  contains  constant  functions, 
and  (6.7),  (6.8),  and  (6.13)  hold. 

(b)  For  all  x  G  such  that  \rif.r\dfl\  >  0,  functions  v  G  Vf  satisfy  = 

0,  and  there  is  a  constant  C,  independent  of  pi  and  v,  such  that 

IklUsfenn)  <  C  (diam(ri’f))  |1  Vw||i,(^.no),  (6.21) 
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for  all  V  €  n  fl)  satisfying  v  =  0  on  dfl.  Moreover  (6.7)  and  (6.8) 

hold  for  u  satisfying  =  0. 

Then  there  exists  iff.  G  Vf  so  that  the  corresponding  function, 


Uap  =  E  G 

satisfies 

1^2 

11^0  -  Wopllfficn)  <  C'(  X!  >  (6-22) 

where  C  is  independent  of  uq  and  v. 

Proof.  We  first  divide  the  set  into  two  disjoint  sets,  namely, 

=  {x  G  Aq  :  rjf.  n  dil  =  0},  and 
A(i'b  =  {xGAa-.v^Pdn^O}. 

Let  iff  G  Vfi,  X  G  Aq,  satisfy  (6.7)  and  (6.8).  Define  (ff,,  for  x  G  Aq /, 
as  in  the  proof  of  Theorem  6.2.  We  know  from  assumption  (a)  that,  for  x  G 
Aq  j,  (6.13)  holds  and  Vfi  contains  constant  functions.  Therefore  following  the 
argument  leading  to  (6.16)  in  Theorem  6.2,  we  get 

Iko  -  ^3LllL2(cna)  <  C'diam(r7^)  e2{x),  x  G  Af^j.  (6.23) 

For  X  G  Aq  we  set  ipf  =  iff.  Now,  =  0  and  from  assumption  (b), 

we  know  that  =  0  for  s  G  Aq  ^.  Thus,  using  (6.21),  with  v  =  UQ  —  iff, 

and  (6.8),  we  have 

IWo  ~  if x\\L2{'n^r]^)  ~  \W  ~ '4’x\\L2{r]^r)^) 

<  C'diam(?7^)e2(x),  xGAf^g.  (6.24) 

Following  the  same  steps  leading  to  (6.17)  in  the  proof  of  Theorem  6.2,  and 
using  (6.23)  and  (6.24),  we  get 

W^o  —  UapW  B2{n)  —  ^  ~  Ilislr/gnO) 

=  C'  X!  Il'^o  “ '^^llL(cno)  +  C'  ^  l!“0  -  ^sllL2(^vnn) 
<  C  (diam(??^))2e2(x).  (6.25) 

xeA^ 

Similarly,  following  the  steps  leading  to  (6.20)  in  the  proof  of  Theorem  6.2,  we 
get 

\\y{u-Uap)\\Yn)<C  Y  ^2^),  (6.26) 
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and  combining  this  with  (6.25),  we  get  (6.22),  where  we  used  the  assumption 
that  diam(?7^)  <  C  for  all  x  €  X''  and  for  all  □ 


Remark  6.2  It  is  clear  from  (6.14)  and  (2.8)  that  if  Uq  is  the  solution  of  (2.1), 
(2.2),  then 


IIUO  -  UGFEM\\H^n) 


<  C'(  ^2(31))  '  , 

xGX‘' 


provided  the  local  approximation  spaces  contain  constant  functions,  and 
(6.13)  holds.  The  above  estimate  is  also  true  if  uq  is  the  solution  of  (2.1),  (2.3) 
provided  conditions  (a)  and  (b)  of  Theorem  6.3  are  satisfied.  We  note  that  in  the 
later  case,  i.e.,  when  uq  satisfies  the  Dirichlet  boundary  condition,  Mo|an  =  0, 
the  space  V^,  corresponding  to  a  particle  x  such  that  r]^  intersects  dfl,  does  not 
need  to  include  constant  functions,  but  the  functions  in  have  to  satisfy  the 
Dirichlet  boundary  condition  on  H  dfl. 


Remark  6.3  The  conditions  (a),  (b)  in  Theorem  6.3,  and  (6.13)  are  known  as 
the  uniform  Poincare  property.  These  conditions  put  restrictions  on  the  shapes 
of  the  For  a  detailed  discussion  on  this  property,  see  [17]. 


Remark  6.4  The  constant  C2  in  (6.3)  is  related  to  the  ratio  of  the  radius  of 
the  largest  ball  contained  in  rjf.  to  the  radius  of  the  smallest  ball  that  contains 
rjf..  A  similar  condition  is  also  assumed  in  the  classical  FEM.  If  this  ratio  is 
uniformly  bounded  for  all  x  €  Aq,  for  all  then  (6.13)  holds. 


Remark  6.5  In  practical  computations,  one  can  easily  construct  particle-shape 
function  systems  (with  k  =  0),  such  that  conditions  (6.2),  (6.3),  (6.13),  and 
conditions  (a),  (b)  of  Theorem  6.3  are  satisfied. 


Remark  6.6  We  observed  that  a  partition  of  unity  is  the  staring  point  for 
the  construction  of  approximating  space  for  GFEM.  It  is  important  to  empha¬ 
size  that  construction  of  partition  unity  for  k  =  0  is  simple,  e.g.,  it  could  be 
constructed  by  Shepard’s  approach  as  discussed  in  Section  4. 


Remark  6.7  We  have  assumed  that  our  particle-shape  function  system  satis¬ 
fies  A1-A7  with  k  =  0  {i.e.,  it  reproduces  polynomials  of  degree  0),  and  have 
seen  that  the  quality  of  the  approximation  in  Theorems  6. 1-6.3  depends  entirely 
on  the  approximability  properties  of  the  spaces  as  quantified  by  ei(x)  and 
€2(3:).  If  we  used  a  particle-shape  function  system  that  reproduced  polynomials 
of  degree  1  {k  =  1),  then  the  space  V'  defined  in  6.6  would  be  enlarged,  and 
its  approximability  would  be  improved,  possibility  only  marginally,  but  this  im¬ 
provement  would  not  be  directly  visible  from  (6.11)  (or  (6.14)  or  (6.18)).  Note 
that  Theorems  6. 1-6.3  are  directed  toward  the  use  of  nonpolynomial  approxi¬ 
mating  functions,  where  the  rate  of  convergence  cannot  be  easily  defined. 

To  clarify  this  point,  suppose  for  (ff  we  use  the  usual  FE  hat  functions  of 
degree  1,  and  Vf  is  the  space  of  constants.  Then  the  GFEM  is  the  classical 
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FEM,  with  the  usual  rate  of  convergence  of  0{h).  However,  (6.11)  (or  (6.14)  or 
(6.18))  does  not  establish  this  rate.  As  a  second  example,  let  be  the  space 
of  linear  polynomials.  Then  the  GFEM  is  a  FE  method,  but  not  a  usual  one. 
The  method  has  the  rate  of  convergence  0(/i^),  but  (6.11)  (or  (6.14)  or  (6.18)) 
only  establishes  0{h). 

More  generally,  if  reproduces  polynomials  of  degree  r,  then  the  func¬ 
tions  which  are  used  in  V‘',  reproduce  polynomials  on  degree  k  +  r.  This 

observation  allows  one  too  establish  the  higher  rate  of  convergence  noted  in  the 
previous  paragraph.  This  will  be  done  in  a  forthcoming  paper. 

The  estimates  in  Theorems  6.2  and  6.3  are  quite  general,  and  allow  us  to 
employ  available  information  on  the  approximated  function  u.  Convergence  of 
the  approximation  can  be  obtained  by  considering  i^i  €  N,  i  =  1,2, . . . ,  such 
that  h'^'  \  0,  where  h''  is  defined  in  (3.80).  This  is  reminiscent  of  the  h- version 
of  FEM.  Convergence  of  the  approximation  can  also  be  attained  by  keeping  v 
fixed,  and  selecting  a  sequence  of  spaces  i  =  1,2,...,  so  that  they  are 

complete  in  H^{fix)  or  in  a  space  V^ifix)  C  that  is  known  to  include  the 

approximated  function  uq.  This  is  a  generalization  of  the  p- version  of  FEM. 

6.2  Selection  of  Vx  and  “Handbook”  Problems 

We  saw  in  Section  6.1  that  it  is  important  to  select  spaces  with  good  local 
approximation  properties.  Principles  for  selecting  shape  functions  that  take  ad¬ 
vantage  of  available  information  on  the  approximated  function  were  formulated 
in  [14,  12].  We  will  use  these  ideas  to  discuss  the  selection  of  the  space  ■  In 
this  section  we  will  suppress  v  in  our  notation. 

Let  iLi(%)  and  iL2(%)  be  two  Hilbert  spaces,  and  suppose  iL2(%)  C 
Then 

dn{H2,Hi)  =  inf  sup  inf  \\u-x\\hi 

S„CHi  U&H2 
dimS„=n  \\u\\„^<l 

in) 

is  called  the  n- width  of  iL2-unit  ball  in  Hi.  Let  Vy  be  an  n-dimensional 
subspace  of  Hi,  and  let 

4'(Pi”\H2,Hi)  =  sup  inf  ||m-x||//,. 
u£H2 
||tl||H2<l 


^{V^\H2,Hi)  is  called  the  sup-inf.  We  will  write  'I'(V(i”^)  for  4'(V]i”\  H2,  Hi) 
if  there  is  no  confusion  about  the  spaces  Hi  and  H2.  It  is  clear  that 

d„(H2,Hi)=  inf  4'(Ei"),H2,Hi). 

dim  —n 

If  an  n-dimensional  subspace  satisfies 

^(Vi”\H2,Hi)  <  C'd„(H2,Hi), 
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where  (7  >  1  is  a  constant,  independent  of  n,  then  we  will  refer  to  as 

a  nearly  optimal  subspace  relative  to  Hi  and  H2.  An  n-dimensional  subspace 
that  satisfies 

^CV^^\H2,Hi)  =  dn{H2,Hi), 

is  referred  to  as  optimal  subspace  relative  to  Hi  and  iJ2-  An  optimal  subspace 
leads  to  the  minimal  error  that  can  be  achieved  with  an  n-dimensional 
space,  namely,  dn{H2,  Hi);  a  nearly  optimal  subspace  leads  to  essentially  the 
same  error,  d„{H2,  Hi). 

Suppose  we  are  interested  in  using  the  GFEM  to  approximate  the  solution 
uo  of  the  Dirichlet  problem, 

J  Auq  =  0  in  n, 

I  uo  =  g,  on  dn, 

where  is  a  bounded  domain  in  Then,  for  each  x  G  X‘',  we  seek  a  finite 
dimensional  space  Vx  that  contains  a  good  approximation  tpx  to  uq  on  rjx  (c/. 
(6.7),  (6.8)).  This  will  be  done  by  taking  advantage  of  the  available  information 
on  uo|^  namely  that  uo\^  is  harmonic.  We  now  illustrate  this  procedure. 

We  now  suppose  that  r/x  is  a  disk  in  and,  for  the  sake  of  simplicity,  suppose 
rjx  is  the  unit  disk.  Let  Hi  =  {u  G  H^^fjx)  '■  u  is  harmonic  in  rjx}.  For  the  space 
H2,  we  use  W{fix),  a  (regularity)  space  known  to  contain  uq.  More  precisely, 
we  suppose  yV{fix)  is  a  linear  manifold  in  {u  G  H^^gx)  ■  u  is  harmonic}  and 
suppose  |||m|||  is  a  norm  on  W{fix)  that  is  rotationally  invariant  and  satisfies 
Il“llffbi7x)  —  II  1^11 1>  “  €  kV’(%).  Moreover,  we  assume  W{f]^  is  complete 

with  respect  to  |||  •  |||,  i.e.,  {W{fi^,  |||  •  |||}  is  a  Hilbert  space.  We  note  that 
yV{Vx)  could  be  any  higher  order  (isotropic)  Sobolev  space. 

It  is  well  known  that  any  u  G  Hi  is  characterized  by  its  trace  on  the  boundary 
I  =  drjx]  these  traces  will  be  in 

S  =  {u{9)  :  0  <  9  <  2tt,u  is  2tt  periodic,  u  G  H^^^{I)}. 

Any  u  G  S  can  be  expanded  in  its  Fourier  series 

00 

u{9)  =  oo  +  cos  k9  +  bk  sin  k9).  (6.27) 

fc=i 

It  is  immediate  that 

00 

l^lffi/2(/)  =  Oq  +  '^i^k  + 

k=l 

where  |u|jji/2(7)  is  a  Sobolev  norm  of  order  1/2  on  I.  So  we  have  a  one-to-one 
correspondence  between  u{r,  9)  G  Hi  ((r,  9)  are  polar  coordinates)  and  u{9)  G  S, 
which  we  express  by  writing  u(r,9)  ~  u{9).  We  easily  find  that 

00 

ll^llffb^x)  “  l''^lffi/2(/)  =  +  ^7)7-  (6.28) 

i=i 
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Thus  we  identify  the  space  Hi  with 

Since  HI'ulH  is  rotationally  invariant,  the  corresponding  norm  on  u{6)  will  be 
translation  invariant,  and  we  thus  can  show  that 

OO 

II|m||P  =  al  +  (6-29) 

j=i 

where,  since  <  lll'^^lll)  we  have  /3j  >  1.  If  we  now  define 

H^{I)  =  {u  G  S  :  \u\i3  <  oo}, 


where 

OO 

Hi  =  al  +  +  bl)k(ik,  (6.30) 

k=l 

then  we  see  that  u{r,6)  e  H2  if  and  only  if  u{9)  e  H^{I)  and  |||m|||  =  \u\13.  We 
thus  identify  the  space  H2  with 

We  will  now  find  an  optimal  subspace  relative  to  Hi  and  H2.  We  will 

exploit  the  correspondence  u{r,  9)  ~  u{9),  and  find  by  first  identifying  an 

optimal  subspace  relative  to  Hi  =  H^^‘^{I)  and  H2  = 

Let  M„  =  {mi,TO2, . . .  ,to„}  be  a  set  of  n  positive  integers,  and  consider 

yM„  _  g  :  u  =  qq  +  ^  (ofc  COS  k9  +  bk  sin  k9)}.  (6.31) 

fceM„ 

Clearly,  is  a  (2n  +  l)-dimensional  space. 

Lemma  6.1  Let  Hi  =  H^^‘^{I),  H2  =  H^{I),  where /3  =  {Pi,  f32,‘ ■  ■),  /3fe  >  1, 
and  let  V^"-  be  as  defined  in  (6.31).  Then 

'l>{V^yH2,Hi)  =  {^{V^-))-"T  (6.32) 

where 

y(y“")  =  inf  A. 

i^M„ 

Proof.  Consider  u  G  H2  given  by 

OO 

u  =  Oo  +  ^^(ofc  cos  k9  +  bk  sin  k9). 

k=l 


Then  from  (6.31),  we  get 


inf  \u 


x\h 


Hi  +  bl)k, 

kGN-Mn 
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where  N  is  the  set  of  all  positive  integers.  Therefore  from  (6.30)  and  the  defi¬ 
nition  of  7(y^"),  we  have 


inf 

xev“" 


_  J2k&N-M„ 

al  +  E  kGN  (al  +  bl)kj3k 

^  J2keN-M„(^k  +  ^k)^ 
J2k&N-M„^^k  +  ^k)^l^k 
1 


Thus, 


sup  inf 
xev«- 


'-X\% 


< 


1 


Let  e  >  0  be  arbitrary.  Then  there  is  an  mo  ^  M„,  mo  >  1,  such  that 


(6.33) 


/3mo  <7(^“")+e- 


(6.34) 


Consider  Umg  =  cosnioO.  Clearly,  m^o  ^  and  therefore  from  (6.27), 


X\hi  ~  Iffi  “  krio. 


Also,  from  (6.30),  we  have  lumol^^  “  ^oPmo-  Therefore,  using  (6.34),  we  get 

^  \Umo-X\s,  1  ^  1 

„6iJ2xev«-  |u|7^  xev«-  \umo\%^  Pmo  x{V^’-)  +  e 
From  this  estimate  and  (6.33),  we  have 


1  .  I  A 1^  i 

7(C“")  +  e  -  xe^“"  “  7(C“") ' 

Since  e  is  arbitrary,  we  get  (6.32).  □ 


Lemma  6.2  Let  Hi  =  and  H2  =  H^{I),  where  f3  =  (/3i,  /92)  •  •  • ))  /3fe  > 

1.  Then 

d2„(i^2,i?i)  =  (7:)-A 

where 

l*n  =  sup  inf  /3,. 

The  proof  of  this  theorem  follows  immediately  from  Lemma  6.1. 
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Theorem  6.4  Suppose  Hi  =  {u  €  H^{'ijx)  ■  u  is  harmonic}  and  H2  =  yV{fix) 
with  the  norm  |||u|||/3  =  \u\p,  given  in  (6.30),  with  I3j>\.  Suppose  in  addition 
that  the  sequence  (ij  is  non-decreasing.  Then  the  space 

=  spanjr-’ cos  j6*,  H  sin  j0}"^Q 

i.e.,  the  span  of  first  (2n  +  1)  harmonic  polynomials  is  optimal  relative  to  Hi 
and  any  H2  (i.e.,  any  of  the  spaces  H2  we  are  considering). 

Proof.  Using  the  correspondence  u{r,  9)  ~  u{9),  we  can  study  the  optimality 
of  a  finite  dimensional  subspace  relative  to  Hi  and  H2,  by  studying  the  opti¬ 
mality  of  a  subspace  relative  to  Hi  and  H2.  The  result  follows  directly  from 
Lemma  6.2.  □ 


Remark  6.8  Obviously  the  condition  on  (3  in  Theorem  6.4  holds  for  any  (iotropic) 
Sobolev  space. 

Remark  6.9  Let  us  return  to  the  solution  of  the  Dirichelet  problem  mentioned 
above.  Suppose  rjx  is  far  from  the  boundary  of  O.  Then  on  rjx,  the  character 
of  the  solution  uq  is  approximately  the  same  in  any  direction.  Thus  it  is  ap¬ 
propriate  to  embed  uq  in  a  space  with  a  rotationally  invariant  norm — a  usual 
(isotropic)  Sobolev  space,  e.g.  And  we  have  learned  that  on  rjx,  uq  is  well 
approximated  by  harmonic  polynomials.  The  situation  is,  however,  somewhat 
different  when  rjx  is  near  the  bounday.  Then  uq  would  be  strongly  influenced  by 
the  boundary  values  g{x).  Hence  some  other  shape  functions,  constructed,  e.g., 
by  the  Handbook  approach,  which  themselves  reflected  these  boundary  values, 
would  be  “best” . 

Thus  the  optimal  shape  functions  are  the  solution  of  the  Laplace  equation. 
This  approach  could  be  also  be  used  for  other  differential  equation,  e.g.,  —Au-\- 
ku  =  0,  or  when  773,  H  U  =  Bi  —  B1/2  where  Bp  is  the  ball  of  diameter  p  and 
homogeneous  normal  boundary  conditions  are  prescribed  on  dTl. 

In  this  section,  we  saw  an  example  of  choosing  an  optimal  local  approxi¬ 
mating  space  Vx,  which  turned  out  to  be  the  span  of  first  (2n  -I-  1)  harmonic 
polynomials.  In  other  problems,  different  local  approximating  spaces,  consisting 
of  optimal  or  near  optimal  approximating  functions,  are  recommended.  These 
optimal  or  near  optimal  approximating  functions  are  solutions  of  other  bound¬ 
ary  value  problems  (posed  on  77^0  U).  Such  locally  posed  problems  are  called 
Handbook  Problems,  and  their  solutions,  which  may  be  available  analytically 
or  computed  numerically,  are  called  Handbook  Functions.  This  nomenclature 
is  reminiscent  of  the  solved  problems  and  their  solutions  (via  formulae,  tables 
etc.)  which  are  used  in  engineering  ([86]).  This  idea  is  also  used  in  commercial 
codes  ([86,  85]). 

One  of  the  main  advantages  of  GFEM  is  that  only  simple  meshes  are  used, 
which  need  not  reflect  the  boundary,  e.g.,  uniform  finite  element  meshes.  Also, 
in  each  rjx,  one  can  use  a  space  Vx  of  arbitrary  dimension  (depending  on  x).  Vx 


could  be  space  of  polynomials  or  any  other  space  of  functions  depending  on  the 
local  properties  of  the  approximated  function. 

Choosing  14  to  be  the  space  of  polynomials  of  low  degree  p  (and  using  {(j)^ 
that  are  reproducing  of  order  k),  we  obtain  the  /i- version  of  FEM.  All  other 
classical  versions  of  FEM — the  p  and  h-p  versions —  are  special  cases  of  GFEM. 
Some  examples  of  the  use  of  this  method  will  be  discussed  in  Section  9. 

7  Solutions  of  Elliptic  Boundary  Value  Problem 

In  this  section,  we  will  discuss  the  approximate  solution  of  the  model  problem 
(2.1)-(2.2)  (or  (2.3)),  introduced  in  Section  2,  by  a  meshless  method.  We  will 
address  the  the  Neumann  boundary  condition  (2.2)  and  the  Dirichlet  boundary 
condition  (2.3)  separately.  These  problems  have  the  variational  formulation 
(2.4). 

For  0  <  /i  <  1,  we  consider  a  family  of  particle-shape  function  systems 

{M^}o<h<i  =  <p^}xex>'}o<h<i, 

satisfying  the  assumptions  A1-A7  in  Section  3.3  and  (3.64).  Recall  that  (3.64) 
is  trivially  satisfied  if  the  shape  functions  are  reproducing  of  order  k.  The  family 
{AI^}o<h<i  was  introduced  in  Section  3.3;  recall  that 

sup  <  h. 

xGX’'  ~ 

We  will  be  interested  in  assessing  the  approximation  error  as  h\0. 

Let  Mo  be  the  solution  of  (2.4),  where  is  a  bounded  domain  with  Lip- 
schitz  continuous  boundary.  Sometimes  in  this  section,  we  will  assume  that 
the  boundary  of  is  smooth.  We  will  use  the  space  defined  in  (3.86), 

to  approximate  uq-  It  was  shown  in  Theorem  3. 1 1  that  is  {k  +  l,q)  reg¬ 
ular.  Moreover,  'Vq\  satisfies  the  local  assumption  LA.  Recall  that  k  is  the 
order  of  the  quasi-reproducing  shape  functions  considered  in  {A4^}  and  q  is  the 
smoothness  index  of  these  shape  functions.  The  parameters  k  and  q  are  in  the 
assumptions  A1-A7  and  we  assume  that  q  <  k  +  1.  We  also  recall  that  'Vq\ 
does  not  involve  all  the  particles  in  X^.  It  only  involves  particles  in  the  set 

=  (7.1) 

Various  classes  of  shape  functions  can  be  used  for  <j)^  in  the  system  {M^}. 
In  Section  4,  one  such  class  of  shape  functions,  namely  RKP  shape  functions, 
were  discussed,  and  references  related  to  other  classes  of  shape  functions  used 
in  practice  were  provided. 

We  note  that  it  is  possible  to  construct  particle-shape  function  system 
satisfying  A1-A7,  such  that  the  set  of  particles  C  fl,  and  the  corresponding 
have  desired  approximation  properties.  We  do  not  consider  such  Vq"),  in 
this  section,  and  we  will  further  remark  on  this  issue  in  the  next  subsection. 
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Let  Us  =  Uh  be  the  approximate  solution  defined  by  (2.7)  with 

S  =  Since  is  {k  +  1,  g)-regular,  we  note  that  'Vq\  C  H  = 

provided  g  >  1.  Thus  Uh  is  the  solution  of 


Uh  G 


Tk,q 

'  a,h 


B{uo,v)  =  /  fvdx,  for  all  u  G 


Tk,q 

'n,h’ 


(7.2) 


where  the  bilinear  form  B  is  either  B,  given  in  (2.5),  or  a  perturbation  of  B. 
Clearly,  Uh  is  the  solution  of  a  Galerkin  method.  This  Galerkin  method  is 
a  meshless  method  since  the  construction  of  the  test  and  the  trial  space,  i.e., 
'Vq\,  does  not  require  a  mesh.  As  we  remarked  in  Section  1,  avoiding  mesh 
generation  is  one  of  the  main  features  and  advantages  of  meshless  methods. 

In  this  section,  we  will  consider  Uh  as  an  approximation  of  Uq  and  primarily 
study  the  error  uq  —  uh-  We  set  some  notations  that  will  be  used  in  this  study 
in  the  following  sections.  We  define 

=  xgaO,  (7.3) 

and 

=  {a-  G  An  :  yf  0}.  (7.4) 

Thus  is  the  set  of  particles  {s}  such  that  f)^  has  non-empty  intersection 
with  on. 


7.1  A  Meshless  Method  for  Neumann  Boundary  Condi¬ 
tion 

In  this  section,  we  will  address  the  approximation  of  solution  uq  of  (2.1)  and 
(2.2)  by  the  meshless  method.  The  analysis  presented  here  is  based  on  the  ideas 
and  results  in  [5,  11].  Also  see  the  references  listed  in  these  articles. 

The  solution  uq  of  (2.1),  (2.2)  can  be  variationally  characterized  by  (2.4), 
which  is 


We  wish  to  approximate  uq  by  Uh,  the  solution  of  (7.2)  with  B  =  B.  For  an 
error  estimate,  from  (2.8)  we  have 

\\uo-Uh\\H^n)<  inf  ||uo  -  xllr/qo)- 
xeK:l 

Suppose  Mo  G  Then,  since  is  (k  +  1,  g)-regular  and  g  >  1,  we  have 

11^0  —  u/illi/qn)  <  Cft,'^||Mo||//!(o), 

where  /x  =  min  (A:,  /  —  1).  We  summarize  this  in  the  following  theorem. 
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Theorem  7.1  Suppose  uq  G  with  I  >  1,  is  the  solution  of  (7.5),  where 

dfl  is  Lipschitz.  Let  Uh  G  with  q  >  I,  be  the  approximate  solution  given 

by  (7.2)  with  B  =  B.  Then 

||uo  -  <  h^\\uo\\H‘{n),  (7-6) 

where 

fj,  =  imn{k,l  —  1)  (7.7) 

We  note  that  the  computation  of  Uh,  in  Theorem  7.1,  depends  on  the  defi¬ 
nition  of  and  involves  particles  that  are  also  outside  LI.  In  the  literature, 
especially  in  the  engineering  literature,  (t,  A:*)-regular  particle  spaces  are  con¬ 
structed  using  particles  inside  LI,  but  the  support  of  some  of  the  corresponding 
particle  shape  functions  could  be  partly  outside  LI.  The  apparent  reason  for 
such  construction  is  that  the  approximate  solution  is  viewed  as  an  interpolant 
with  respect  to  data  inside  LI,  and  hence  the  particles  that  are  only  inside  LI  are 
considered.  This  certainly  is  not  necessary. 

The  construction  of  the  approximation  space  S  (in  (2.7))  using  particles 
only  inside  LI  sometimes  may  lead  to  better  conditioning  of  the  underlying  lin¬ 
ear  system.  On  the  other  hand,  such  construction  is  more  expensive  and  the 
approximations  could  show  boundary  layer  behavior  ([13]). 

7.2  Meshless  Methods  for  Dirichlet  Boundary  Condition 

In  this  section,  we  consider  the  approximation  of  the  solution  uq  of  the  Dirichlet 
boundary  value  problem  (2.1)  and  (2.3)  by  meshless  methods.  The  variational 
characterization  of  uq  is  given  by 

j  uoG  Hq(LI) 

B(uq,v)  =  J  fvdx,  for  all  u  €  i7g(n).  ^  ^ 

The  Galerkin  method  (2.7)  to  approximate  uq  would  require  that  the  approxi¬ 
mating  space  S'  be  a  subspace  oi  H  =  Hl(Lt)  and  thus  that  the  approximating 
functions  satisfy  the  essential  homogeneous  Dirichlet  boundary  condition.  Un¬ 
like  shape  functions  used  in  FEM,  the  particle  shape  functions  (jff  (we  consider 
h  as  the  parameter),  considered  in  Section  3.3,  do  not  in  general  satisfy  the  so 
called  “Kronecker  delta”  property,  i.e.,  (j)^(y)  yf  bx,y,  x,y  G  X^.  This  is  also 
true  for  translation  invariant  particle  shape  functions  discussed  in  Section  3.2 
(see  Section  4.2).  Thus  it  is  difficult  to  construct  a  subspace  S  C  such 
that  S  could  be  used  in  (2.7)  as  the  approximation  space  and  the  functions  in 
S  satisfy  the  Dirichlet  boundary  condition. 

In  the  literature,  several  meshless  methods  have  been  proposed  to  approx¬ 
imate  the  solutions  of  Dirichlet  boundary  value  problems.  They  are  meshless 
methods  in  the  sense  that  they  use  as  the  approximating  space.  These 
methods  are: 


91 


1.  The  Penalty  Method 

2.  The  Lagrange  Multiplier  Method 

3.  The  Nitsche  and  Related  Methods 

4.  The  Collocation  Method 

5.  Combination  of  meshless  and  finite  element  method 

6.  The  characteristic  function  method 

In  this  section  we  will  describe  these  methods.  We  note  that  GFEM,  dis¬ 
cussed  in  Section  6,  uses  an  approximating  space,  different  from  'Vq\,  and  can 
also  be  used  to  approximate  the  solution  of  a  Dirichlet  boundary  value  problem. 

We  will  assume  that  the  boundary  d^l  of  is  sufficiently  smooth.  The 
smoothness  assumption  on  the  boundary  simplifies  the  arguments  presented 
here,  but  various  results  could  be  obtained  when  the  boundary  is  not  smooth. 


1.  The  Penalty  Method. 

The  main  idea  of  the  penalty  method  is  to  use  a  perturbed  variational  prin¬ 
ciple.  For  CT  >  0,  we  consider  the  bilinear  form 

B{u,  v)  =  B„{u,  v)  =  B{u,  v)  +  h~'^D{u,  v),  (7.9) 

where 

B{u,v)  =  /  {Wu  ■  Wv  +  uv)  dx,  (7.10) 

Jn 

D{u,v)  =  /  uvdx.  (7.11) 

Jdn 

We  note  that  (7.10)  is  the  bilinear  form  given  in  (2.5).  We  consider  the  solution 
Uh  =  Ucr,h  G  of  (7.2),  namely 


Ba{ucr,h,v)=  [  fvdx,  for  all  V  G  (7.12) 

Jn 

We  note  that  is  us,  where  us  is  defined  in  (2.7).  For  v  G  let 

Qcr{v)  =  B{v,v)  +  h~'^ D{v,v)  —  2  f  fvdx.  (7.13) 

Jn 


It  is  well  known  that 


Qa{ua,h)  =  min  Qa{v). 

■■ 


^  Q,,h 


(7.14) 


We  now  present  a  convergence  result  for  the  penalty  method. 


Theorem  7.2  Suppose  Ug  G  B‘(fl)  n  BQ(fl),  I  >  3/2,  is  the  solution  of  (7.8). 
Let  Ua,h  G  be  the  solution  of  (7.12).  Then  for  any  0  <  e  <  min(l  — 

3/2,  1/2),  we  have 

lluo  —  Uo-./i||ffi(n)  <  C{e)h^\\ug\\Hi(a),  (7-15) 
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where 

^  =  min(fc,  ;  -  1,  -,  /c  +  -  -  -  -  e,  /  -  -  -  -  -  e),  (7 

and  C(e)  depends  on  e,  but  is  independent  of  h  and  uq. 

Proof.  For  any  v  G  we  define 

Rcr{v)  =  B{uo  —  v,uo  —  v)  +  h~'^ +  v,  +  v).  (7 

on  on 

Then  from  Green’s  Theorem, 

Ra{v)  =  B{uo,uo)  +  B{v,v)  -  2B{uo,v) 

+B{v,v)  +  h~'^ D{v,v)  —  2  /  fvdx 

Jn 

=  B{uo,uo)  +  +  Qcr{v),  for  all  G 

where  Qa{v)  is  given  by  (7.13).  Therefore, 

min  R^iy)  =  B{uQ,UQ)  +  h'^D{^^,^^)+  min  Qa{v), 

-ev^a’X  -evkl 


and  thus  from  (7.14)  we  get. 


Ra{ucr,h)  =  min  Rcy{v) 


Hence,  from  (7.17)  and  the  above  relation,  we  have 

11^0  II //i (n)  B{tiQ  Ufj  fi^  uq  Ufj  fi) 

^  Rfjip^a.h) 

<  Ra{v),  for  all  G  .  (7.18) 

Since  is  (fc  +  1,  (7)-regular  with  g  >  1,  there  is  gt  G  such  that 

llwo  -  gh\\Hyn)  <  C'^^||uo||H*(n),  (7-19) 

where  g  =  min(fc  +  1  —  s,l  —  s)  and  0  <  s  <  1.  Now  from  (7.17)  with  v  =  gh 
and  using  Cauchy-Schwartz  inequality  we  have 

Ra{gh)  <  C  ^||mo  -  ghWnyn)  +  h'^ 
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We  will  estimate  the  right  hand  side  of  the  above  inequality.  We  first  note  that 
t6o  =  0  on  (9n.  Let  0  <  e  <  min(Z  —  |,  ^).  Then  using  a  trace  inequality  and 
(7.19)  with  s  =  (1/2)  +  e,  we  get 

Wgiillt^idn)  =  W'^o  —  9h\\L2(da)  ^  C'(e)ll'Wo  — 

<  C{e)h^^^\\uorH.^a),  (7-21) 

where  /xi  =  min(fc  +  ^—  e,  I— ^  —  e).  Also  from  a  trace  inequality,  we  have 

ll^lli.,»,<C(c)||«o||S„^„,.  (7.22) 

Now  using  (7.21),  (7.22),  and  (7.19),  with  s  =  1,  in  (7.20),  we  get 

RA9k)  <  C{e)  ^  ^  Iholl^^qn) 

<  C(e)/i2A^||Mo||^q^),  (7.23) 

where  ^  =  min(fc,  I  —  1,  |,  fc+  ^  |  —  e,  I  —  |  |  —  e).  Finally,  combining  (7.18) 

and  (7.23)  we  get  the  desired  result.  □ 

Remark  7.1  If  we  consider  in  Theorem  7.2  such  that  k  +  I  >  I  >  3/2, 
then  with  cr  =  1  ^  —  e,  it  is  easy  to  see  that  (7.15)  holds  with 

.  =  ^a-^-e). 

The  estimate  (7.15)  can  be  improved.  We  present  the  following  result,  based 
on  the  analysis  in  [4,  5,  6,  8],  without  proof. 

Theorem  7.3  Suppose  Uq  €  i7*(n)  n  TLq  (fl)  is  the  solution  of  (7.8).  LetUa^h  G 
Vq’^  be  the  solution  of  (7.12).  If  k  +  1  >  I  >2,  then,  for  any  e  >  0,  we  have 

||mo  -  Ua,h\\H^(o.)  <  C{e)h^~''\\uo\\Hi(n),  (7.24) 

where  C{e)  is  independent  of  Uq  and  h  but  depends  on  e,  and  p,  is  given  by 

p  =  min  ^cr,  ;  +  cr  -  2,  ;  +  I  -  _  i)^  ^  (7,25) 

where 

K  =  max(l,  — (,  (7.26) 

Remark  7.2  If  cr  >  1  in  Theorem  7.3,  then  l  +  a  —  2  >  /  +|  —  |,  and  therefore 
p  in  (7.24)  is  given  by 

/i  =  min  ^cr,;+ I  -  — -{I  -  1)^  ,  (7.27) 

where  n  is  as  in  (7.26). 
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Example:  Consider  I  =  2  and  ct  =  1  in  Theorem  7.3.  Then  from  (7.26),  we 
have  K  =  1,  and  (7.27)  yields  /i  =  min(l,  1, 1)  =  1.  This  is  the  optimal  rate  of 
convergence.  For  higher  values  of  I  and  k  +  1  >  I,  there  is  a  loss  in  the  rate  of 
convergence  and  we  get  sub-optimal  rate  of  convergence,  i.e.,  /i  <  min(/  —  1,  k). 

Thus  from  Theorem  7.3,  we  conclude  the  following: 

•  It  is  advantageous  to  use  with  higher  values  of  k  since  it  leads  to 
higher  accuracy.  For  example,  if  I  =  4  and  cr  =  3,  then  k  =  2  and  (7.27) 
yields  /i  =  3(^^).  Thus  higher  values  of  k  will  increase  accuracy.  But 
higher  values  of  k  reduce  the  sparsity  of  the  resulting  linear  system. 

•  Too  small  or  too  large  value  of  a  may  decrease  the  accuracy  significantly. 
For  example  if  cr  =  2fc  -|-  1,  then  k  =  k  +  1  and  7.27)  yields  /i  =  0. 

The  use  of  penalty  method  was  recently  suggested  in  the  literatue,  e.g., 
[3],  [57],  [89],  without  any  theoretical  analysis.  A  concrete  penalty  value  of  a, 
unrelated  to  I  or  k,  was  suggested  in  [89]  based  on  experience. 

2.  The  Lagrange  Multiplier  Method 

The  theory  of  Lagrange  multiplier  method,  in  the  context  of  finite  element 
method,  was  developed  in  [7]  (see  also  [11]).  This  theory  can  also  be  extended 
to  meshless  methods. 

It  is  known  (c/.  [11])  that  the  effectiveness  of  this  method  depends  on  a 
delicate  relation  between  the  approximating  space  (12)  and  the  space  of 
Lagrange  multipliers  S'^’*  (512),  where  both  (12)  and  (512)  satisfy  an 
inverse  assumption.  In  the  context  of  meshless  methods,  we  let  the  approxi¬ 
mating  space  to  be  the  particle  space  V^'^.  We  know  that  is  (fc  -I-  !,(?)- 
regular,  and  satisfies  the  inverse  assumption,  lA,  under  additional  hypotheses 
given  at  the  end  of  Section  3.3.  The  space  of  Lagrange  multipliers  5'^['’^’'^(512) 
has  the  same  (2,  A:*)-regularity  as  the  approximating  space,  and  the  functions 
in  5'^['’^’*(512)  are  defined  only  on  512  with  respect  to  particles  on  512.  Thus, 
512  must  contain  enough  particles.  We  note  that  the  functions  in  5'^['’^’^(512) 
are  not  restrictions  of  functions  in  on  512.  Then  following  the  analysis 
in  [7,  11],  one  can  show  that  if  the  size  of  the  supports  of  the  basis  elements 
of  £'^['’^’*(512)  is  of  the  same  order  as  \E^\,  x  €  {E^  and  Aq^  defined 

in  (7.3)  and  (7.4),  respectively),  then  the  approximate  solution  obtained  from 
the  Lagrange  multiplier  method  converges.  If  the  size  of  the  supports  of  the 
basis  elements  of  £^[''^’'^(512)  is  smaller  than  r]^,  x  G  Ag^,  then  the  method  is 
unstable.  This  relationship  was  further  analyzed  in  [72]  and  [73]. 

The  Lagrange  multiplier  technique  leads  to  the  optimal  rate  of  convergence 
in  comparison  to  the  penalty  method,  where  the  optimal  rate  of  convergence  is 
usually  not  attained.  But,  as  we  mentioned  before,  the  sufficient  conditions  for 
convergence  are  quite  delicate. 

Recently,  the  Lagrange  multiplier  technique  was  applied  in  the  context  of 
meshless  methods  without  any  theoretical  analysis  in  [24],  [57],  [60],  [67]. 
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3.  The  Nitsche  and  Related  Methods 

Because  of  the  delicate  nature  of  Lagrange  multiplier  method,  there  was  an 
interest  to  look  for  other  methods  to  deal  with  the  issue  of  the  imposition  of 
Dirichlet  boundary  conditions  and  to  avoid  complications  that  are  present  in 
Lagrange  multiplier  method.  Towards  this  end,  certain  methods  were  proposed 
in  [20]  and  [79].  But  much  earlier,  a  similar  method  was  proposed  by  Nitsche 
in  [70].  We  will  discuss  the  Nitsche  method,  following  the  presentation  in  [79]. 
We  will  still  assume  that  dfl  is  smooth. 

To  approximate  the  solution  uq  of  (7.8)  by  Nitche  method,  we  consider  the 
particle  space  'Vq\  with  q>2.  We  also  assume  that 

•  Card  {Ag^)  <  k  and 

Cih  <  hEH  <  C2h,  X  G  (7.28) 

where  h^h  =  \E!^\,  for  x  G  Aq^^-,  and  A^^  are  defined  in  (7.3)  and 
(7.4),  respectively. 

•  There  is  0  <  /C  <  oo,  /C  =  such  that 

II  1^11- 1, ft  <  /C[S(u,u)]l/^  for  all  v  G  (7.29) 

where  B(u,v)  was  defined  in  (7.10)  and 


We  define  the  bilinear  form  B{u,v)  =  B^(u,v),  where 

dxL  dv  f 

B^{u,v)  =  B{u,v)  —  D{  —  ,v)  —  D{  —  ,u)  +  j  ^  /  uvds, 

UTi  UTl  —  J 


with  7  >  0;  B{u,v),  D{u,v)  are  as  defined  in  (7.10),  (7.11),  respectively.  The 
approximate  solution  Uh,'y  G  obtained  from  the  Nitche  method,  is  given 

by 

B^{uh^-y,v)  =  (  fvdx,  for  all  u  G  (7.31) 

Jn 


For  u  G  77^  (n),  we  define  the  norm 


Fin 

HIP  =  i?(u,u)  +  ||-f_._,+ hill.,. 


where 


IffO(Sh), 
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and  II  11^1  ^  is  given  by  (7.30)  with  v  replaced  by  u.  We  first  note  that  from 


Schwartz  inequality,  we  have 

h'^  f  uvds  < 


1/2 


1/2 


/S'* 


ds 


lE>' 


ds 


< 


316  ^af2 

(7.32) 


=  \\u\\i,h\\y\\ih- 


Also, 


dv 


^  /  C/t'  s  X — ^ 


9n 


<  IklU, Jl/j-ll-i.- 


dn 


l-l 

'an' 


Using  the  same  arguments  used  to  obtain  (7.33),  we  get 


^,du  ,  ,,du,, 


\k,h- 


Now  using  (7.32)-(7.34),  it  is  can  be  easily  shown  that 

Bj{u,  w)  <  (1  +  7)  II  |u||  I  II  liill  I- 


(7.33) 


(7.34) 


(7.35) 


We  now  show  that  for  a  proper  value  of  7,  the  bilinear  form  Bj(u,v)  is 
coercive. 


(7.36) 


Lemma  7.1  Suppose  K.^  <  7,  where  K.  satisfies  (7.29).  Then, 
i?yu^i)>C'r|||u|||2,  forallveY% 
where  C*  =  C*{X^)  >  0. 

Proof.  Let  and  let  e  >  0  arbitrary.  From  the  definition  of  Bj{u,v) 

and  (7.33)  with  u  =  v,  we  have. 

Fin 

Bj{v,v)  =  B{v,v) -2D{v,—)  +  j\\v\\lj^ 


Fin 

>  i?(.,n)-2||n||.,,||-||_.,,+7lM||,, 


>  B(i;,r;)  -ellnliy^  -  ^ II II?. 4 +  bikll i. 


1  Fin 

=  i?(.,n)--||-|y._,  +  (7-e)||u|||_ 

>  (1  -  ^)B(i;,i;)  +  (7- e)||i;|R 
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(7.37) 


Therefore,  considering  e  =  +  7)  in  the  above  inequality,  we  get 

B^{v,v)  >  Ci[B{v,v)  +  ||w|||_^], 

where  Ci  =  min()^~^2  >  ^ — )•  Now,  from  the  definition  of  |||  •  |||  and  using 
(7.29),  we  get 

111^^111^  <  .B(x;,  v)  + /C^i?(v,  v)  +  ||u|||  ^ 

<  {l  +  IC^)[B{v,v)  +  \\v\\lJ. 

Thus  combining  the  above  inequality  with  (7.37)  we  get 

B,{v,v)>C*\\\v\\\^ 

where  C*  =  Yq;^min(^^^,  ),  which  is  (7.36).  □ 

We  now  present  the  following  result: 

Theorem  7.4  Suppose  uq  G  H^Q),  I  >  2  is  the  solution  of  (7.8).  Let  Uh,^  G 
'Vq\,  with  q>  2,  be  the  solution  of  (7.31),  where  satisfies  (7.28),  (7.29). 
Then 

\\uo  -  Uh,y\\m{n)  <  ^  =  min(A:,;  -  1),  (7.38) 

where  C*(X‘')  is  as  in  (7.36). 

Proof.  It  is  easy  to  see  that 

B.y{uo,v)  =  [  fvdx,  for  all  V  G  i7^(fl), 

Ja 


and  therefore, 

B.y{uo  —  Uh,-y,  v)  =  0,  for  all  v  G 
Now  for  any  gh  G  using  (7.36),  (7.39),  and  (7.35)  we  have 

II  ln/i  ^/i,7  III  ^  Bj  i^gh  Ufi^j  ^  9h  ^71,7) 

—  Bj{gfi  Uq,  gii  Ufij) 

<  ^^IIK-5/.||||||5.-z^..7lll, 

and  hence 

(1 

\\\9h  -  uh,j\\\  <  IIIuq  -  gh\\\- 


(7.39) 
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Therefore, 


IIIWO  -  ■U/i,7l||  < 


<  CII |uo  -gh\\\,  for  all  g  € 


Now  using  (7.28)  and  a  trace  inequality,  we  have 

<  Ch  ^\\uo  —  gh\\‘HO(^QQ'f 
~  ^  —  5/illHO(n)  + /iIIuq  — 

=  C  (ji  ^||wo  —  5/illffO(Q)  +  ||wo  —  5h|llf-i(n))  ,  (7.41) 

where  C  depends  on  k.  Also  using  a  similar  argument,  we  have 

II <  C  [h^uo  -  +  ho  -  3.111,1(0))  •  (7.42) 

where  C  depends  on  k.  Thus  from  (7.41)  and  (7.42)  we  get, 

llko  — 3.IIP  =  ho  —  3.111,1(0)  +  II  ^  ^  ^  111.  1^,1  +  ho  —  3/i|||_/i 

2 

i=o 

and  hence  from  the  definition  of  |j|  •  |||  and  (7.40),  we  have 


ho  —  h,7lhi(0)  <  II  1^0  —  11,1,7 

2 


< 


z 

^/i^ho  -  3.lhi(a), 
j=o 


for  all  gh  G  V^(h 


Finally,  we  choose  g  G  'Vq\  such  that 


(7.43) 


\\uo  —  g\\H‘(ii)  ^  *holh‘(0))  0  <  s  <  2, 

where  =  min(fc+  1,1).  Using  this  in  (7.43)  we  get  the  desired  result. □ 

We  now  discuss  the  situations  where  the  assumption  required  to  prove  The¬ 
orem  7.4  are  satisfied.  The  major  problem  is  to  estimate  JC{X^)  given  in  (7.29). 
We  would  like  to  have  JC{X^)  <  C,  uniformly  for  all  0  <  /i  <  1.  If  the  supports 
77(1  of  the  particle  function  (f)’^  are  “reasonable”,  e.g.,  circles  in  or  spheres  in 
or  similar,  then  it  is  easy  to  see  that  the  necessary  condition  for  /C(A^)  <  C 
is  that,  hEh  >  alry^jl  for  x  G  A|q.  This  can  be  enforced  by  properly  selecting  the 
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set  of  particles  This  aspect  can  also  affect  the  design  of  adaptive  meshless 
(Nitsche)  method.  Since  the  estimates  of  these  constants  difficult  to  estimate 
accurately,  we  may  select  larger  value  of  7  in  (7.31)  so  that  the  Theorem  7.4  is 
valid. 

The  Nitsche  method  presented  here  is  superior  to  both  the  penalty  method 
and  Lagrange  multiplier  method.  Nitsche  method,  in  the  framework  of  meshless 
methods,  was  addressed  in  [15]  and  implemented  in  [76]  and  [46]. 

4.  The  Collocation  Method 

Collocation  method,  in  the  framework  of  meshless  methods,  was  recently 
proposed  in  [3,  89,  47].  The  method  consists  of  adding  constraint  equation, 
at  certain  points  of  the  boundary  dfl,  to  the  stiffness  matrix.  No  analysis  was 
presented  to  address  the  convergence  of  the  approximate  solution  obtained  from 
this  method. 


5.  Combination  of  Meshless  Methods  and  the  Finite  Element 
Method 

This  method  was  proposed,  e.g.,  in  [51].  The  main  idea  in  this  approach 
is  to  use  classical  finite  elements  (which  could  also  be  interpreted  as  particle 
functions)  in  a  neighborhood  of  the  boundary  dfl,  and  to  select  other  particle 
functions  such  that  their  supports  do  not  intersect  dfl. 


6.  Characteristic  Function  Method 

The  method  was  proposed  in  connection  to  Ritz  method  when  the  approx¬ 
imating  functions  were  global  polynomials  (see  [63], [50]).  If  a  domain  has  a 
smooth  boundary  dfl,  there  exists  a  smooth  function  4)  such  that 


4)(a;) 

4)(a:) 

and  |V4)(x)|  >  a 


>  0,  a;  G  fl, 

=  0,  a;  G  dn, 

>  0,  a;  G  dn. 


Let  S'^  =  {u  :  u  =  4)u,  v  G  Then  it  is  obvious  that  C 

We  approximate  the  solution  uq  of  (2.1)  and  (2.3)  by  Uh  G  5”^,  where  Uh  is  the 
solution  us  of  (2.7)  with  S'  =  5')f. 

For  uq  G  H^{Vi)  n  i7g(n),  /  >  2,  we  define  wq  =  Then  using  Hardy’s 
inequality  (Thm.  329  of  [49]),  one  can  show  that  wq  G  77^“^  (H).  Using  this 
result,  we  obtain  the  following  theorem. 


Theorem  7.5  Suppose  uq  G  i7*(U)  n  Hq{Q),  and  suppose  I  >  2.  Then  there 
exists  Wfi  G  such  that  gt  =  ^Wh  satisfies 


\\uo-gh\\m{<P)<Ch^\\uo\\H‘{n),  p,  =  -  2).  (7.44) 
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Proof.  Recall  that  'Vq\  is  {k  +  1, (7)-regular  with  g  >  1.  Then  there  exists 
Wh  G  such  that 

ll'ii'o  —  Wh\\H^{n)  <  Ch^\\wQ\\Hi-i{a)  <  C'/i^||'Uo||//!(a),  (7-45) 

where  fx  =  mhi(k,  1  —  2).  Now  from  the  definition  of  Wq,  we  have 
Uo  -  ^Wh  =Uo-  ^Wo  +  ^{Wo  -  Wh)  =  ^{Wo  -  Wh), 
and  hence,  using  (7.45),  we  have 

||mo  -  <  C\\wo  -  M  =  min(A:,/  -  2), 

which  is  the  desired  result.  □ 

Remark  7.3  It  is  clear  from  (2.8)  and  (7.44)  that  for  I  >  2, 

11^0  -  UhWH^n)  <  Ch^\\uo\\H‘(n),  k-  =  mm{k,l-  2). 

We  further  note  that  this  order  of  convergence,  or  (7.44)  in  the  last  theorem, 
cannot,  in  general,  be  improved. 


7.  The  Generalized  Finite  Element  Method 

We  note  that  all  the  methods  described  so  far  primarily  use  as  the 
approximating  space.  The  GFEM,  on  the  other  hand,  uses  different  approxi¬ 
mating  spaces,  as  we  have  seen  in  Section  6.  The  use  of  these  approximating 
spaces  makes  GFEM  extremely  flexible. 

We  recall  that  in  GFEM,  we  start  with  a  partition  of  unity  with  respect 
to  a  simple  mesh  that  need  not  conform  to  the  boundary  of  the  domain.  This 
partition  unity  could  be  the  particle  shape  functions  defined  in  Section  3.  Then 
“handbook”  functions  are  used  as  local  approximating  spaces.  The  Dirichlet 
boundary  condition  can  be  implemented  by  choosing  the  local  approximation 
space  14)  for  x  €  4gQ,  such  that  the  functions  in  14  satisfy  the  Dirichlet  bound¬ 
ary  conditions. 

We  presented  a  few  approaches  on  how  to  use  meshless  approximation  to 
approximate  solutions  of  PDE’s.  To  impose  Dirichlet  boundary  conditions  on 
meshless  approximation  is  a  challenge,  and  we  looked  into  some  methods  that 
can  overcome  this  difficulty.  While  discussing  these  methods,  we  assumed  that 
the  boundary  of  the  domain  is  smooth  for  simplicity.  But  the  results  presented 
here  can  be  generalized  to  include  non-smooth  boundary,  especially,  piecewise 
smooth  boundary. 

Some  of  methods  were  implemented  and  reported  in  the  literature,  but  they 
lacked  rigorous  theoretical  analysis.  All  the  methods  reported  here  have  certain 
advantages  as  well  as  disadvantages.  If  the  particle  space  is  used  as  an 
approximating  space,  in  our  opinion,  the  Nitsche  method  is  very  promising 
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because  it  is  robust  relative  to  other  methods  and  is  easy  to  implement.  But 
we  note  that  is  difficult  to  construct  for  higher  values  of  k,  and  the  use 
of  with  lower  values  of  k  reduces  the  accuracy  of  the  method.  On  the 

other  hand,  GFEM  uses  a  partition  of  unity  with  A:  =  0),  which  is 

easy  to  construct,  and  higher  accuracy  can  be  attained  by  using  suitable  local 
approximation  spaces. 


8  Implementational  Aspects  of  the  Meshless  Method 

In  this  section,  we  will  briefly  discuss  the  implementational  aspects  of  meshless 
methods  and  the  GFEM.  Similar  to  the  finite  element  method,  the  implemen¬ 
tation  of  meshless  methods  and  the  GFEM  has  four  major  parts: 

1.  Gonstruction  of  particle  shape  functions. 

2.  Gonstruction  of  the  stiffness  matrix. 

3.  Solution  of  the  linear  system  of  equations. 

4.  A-posteriori  error  estimation,  adaptivity,  and  computation  of  data  of  in¬ 
terest. 

We  now  discuss  these  items  in  turn. 

1.  Construction  of  particle  shape  functions 

In  the  classical  finite  element  method,  one  starts  with  a  mesh  that  is  related 
to  the  domain,  and  then  shape  functions  are  defined  with  respect  to  the  cho¬ 
sen  mesh.  In  a  meshless  method,  one  starts  with  particles  {x}.  Gorresponding 
to  each  particle  x,  a  particle  shape  function  with  compact  support  is  con¬ 
structed,  such  that  rjaj’s  form  an  open  cover  of  the  domain  The  construction 
of  shape  functions  that  are  reproducing  of  order  fc  =  0  or  1  is  not  difficult. 

For  fc  =  0,  one  may  use  Shepard’s  approach  ([77])  as  described  in  Section  4.1. 

For  k  =  1  and  for  an  appropriate  particle  distribution,  one  may  first  construct 
a  mesh  using  tetrahedrons  such  that  the  particles  are  the  nodes  of  the  mesh 
(i.e.,  the  vertices  of  the  tetrahedrons).  This  procedure  is  not  difficult  as  there 
are  efficient  codes  available  for  constructing  such  a  mesh.  The  shape  function 
corresponding  to  the  particle  x  can  be  taken  to  be  the  standard  hat  functions, 
whose  support  is  the  union  of  all  the  simplices  with  x  as  one  of  its  vertices.  We 
note  however  that,  for  k  =  1,  smoother  shape  functions  can  also  be  constructed 
(see  [48]).  For  fc  =  0,  1,  one  has  to  check  that  card(5'^)  <  k,  «:  is  independent 
of  X,  where  Sx  =  {y  ■  Vy  Vx  0}-  For  the  Nitsche  Method,  described  in 
Section  7.2,  one  also  has  to  check  that  1C,  defined  in  (7.29),  is  bounded.  The 
construction  of  particle  shape  functions  for  fc  >  2  is  more  expensive  than  for 
k  =  0, 1,  and  it  may  be  more  difficult  to  check  assumptions  A1-A7  and  (3.64), 
which  ensure  convergence. 
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In  contrast,  the  GFEM  uses  only  a  partition  of  unity,  and  thus  particle 
shape  functions  with  A:  =  0, 1,  described  in  the  last  paragraph,  can  be  used  for 
this  purpose.  Also,  a  simple  regular  distribution  of  particles  could  be  used  to 
construct  the  partition  of  unity.  The  space  of  local  shape  functions,  V^,  could 
be  created  analytically  or  through  “handbook”  solutions.  Dirichlet  boundary 
condition  is  also  treated  by  an  appropriate  selection  of  14,  and  hence  one  does 
not  have  to  use  special  methods,  e.g.,  the  penalty  method,  Nitsche  method,  etc., 
which  simplifies  the  implementation. 

2.  Construction  of  the  stiffness  matrix 

The  construction  of  the  stiffness  matrix  for  a  meshless  method  is  laborious 
and  delicate.  In  fact,  this  is  where  one  pays  the  price  for  avoiding  mesh  gener¬ 
ation.  The  elements  of  the  resulting  stiffness  matrix  are  integrals,  which  have 
to  be  numerically  evaluated  over  various  regions.  These  regions  are  not  simple 
tetrahedrons  as  in  the  finite  element  method,  where  they  naturally  come  from 
a  mesh.  These  regions,  for  a  meshless  method,  are  of  the  form  rjx  H  rjy  O  fl, 
x,y  G  X‘',  and  can  be  extremely  complicated.  Also  the  integrals  have  to  be 
evaluated  accurately  as  it  is  known  that  inaccurate  numerical  integration  leads 
to  very  poor  results  (see,  for  example,  [28]).  A  special  numerical  integration 
scheme  is  given  in  [30],  where  are  spheres  and  the  region  of  integration 
are  the  intersection  of  two  spheres.  The  problem  of  effective  integration  has 
also  been  addressed  in  [33,  43,  76,  82,  83].  The  numerical  integration  poses 
additional  problems  in  GFEM  when  singular  functions  are  included  in  the  local 
approximating  space  14.  Standard  integration  schemes  in  this  situation  yield 
poor  accuracy.  This  problem  in  GFEM  was  handled  in  [82]  by  using  adaptive 
numerical  integration.  Because  of  this  sensitivity  to  numerical  integration,  the 
use  of  adaptive  integration  is  preferred,  in  general,  in  GFEM. 

Thus  we  see  that  an  accurate  and  effective  numerical  integration  scheme  to 
approximate  the  elements  of  the  stiffness  matrix  is  essential  for  the  success  of 
meshless  methods.  We  will  further  remark  on  this  issue  in  the  next  item  of  this 
section.  We  note  that  numerical  integration  and  construction  of  stiffness  matrix 
in  these  methods  are  parallelizable. 

3.  Solution  of  the  linear  system 

The  exact  stiffness  matrix  (without  numerical  integration)  obtained  from  a 
meshless  method  could  be  positive  definite  with  a  large  condition  number.  This 
is  caused  by  using  shape  functions  with  large  overlap  between  their  supports, 
which  makes  the  shape  functions  “almost”  linearly  dependent.  Moreover,  the 
exact  stiffness  matrix  obtained  from  GFEM  could  be  positive  semi-definite,  as 
shown  in  [82].  But  the  underlying  linear  system  obtained  from  the  GFEM  is 
always  consistent,  i.e.,  the  linear  system  has  non  unique  solutions.  The  lack  of 
unique  solvability  of  the  linear  system  does  not  imply  that  the  GFEM  produces 
non-unique  solutions.  In  fact,  if  the  vector  {cx,j}i<j<n^,  with  dim  14  =  nx,  is 
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a  solution  of  this  linear  system,  then  the  solution, 

Ufi  =  'y '  y '  Cx,j  V’k  j 

X  j  =  l 

obtained  from  the  GFEM,  where  tpl  is  a  basis  of  t4,  is  unique. 

We  already  have  mentioned  the  importance  of  numerical  integration  in  eval¬ 
uating  the  elements  of  the  stiffness  matrix  obtained  from  a  meshless  method. 
We  further  note  that  the  elements  of  the  load  vector  is  also  evaluated  by  numer¬ 
ical  integration.  To  obtain  a  consistent  linear  system  (after  the  use  of  numerical 
integration),  the  numerical  integration  scheme  applied  to  compute  an  element  of 
the  load  vector  should  be  same  as  the  scheme  used  to  compute  the  corresponding 
row  of  the  stiffness  matrix. 

To  find  the  solution  of  the  linear  system  obtained  from  a  meshless  method 
(or  from  the  GFEM) ,  one  can  use  a  special  direct  solver  based  on  elimination  or 
an  iterative  solver.  In  [82],  direct  solvers,  e.g.,  subroutines  MA27  and  MA47  of 
Harwell  Subroutine  Library,  was  used  to  solve  the  sparse  positive  semi-definite 
linear  system  obtained  from  GFEM.  The  use  of  these  solvers  was  successful  even 
when  the  nullity  of  the  stiffness  matrix  was  large.  It  was  also  shown  in  [82]  that 
round-off  errors  did  not  play  a  significant  role  in  solving  the  linear  system,  i.e., 
the  round-off  error  was  almost  same  as  when  the  standard  finite  element  linear 
system  is  solved  by  the  elimination  method. 

An  iterative  algorithm  for  solving  such  linear  systems  was  given  in  [82]. 
The  idea  of  this  algorithm  is  to  perturb  the  stiffness  matrix  by  a  unit  matrix 
multiplied  by  a  small  parameter.  The  perturbed  matrix,  say  P,  is  positive 
definite  and  any  solver  could  be  used  to  solve  Px  =  b.  Using  this  fact  and  a  few 
iterations  of  a  simple  iterative  technique,  a  solution  of  the  original  linear  system 
could  be  obtained.  We  refer  to  [82]  for  a  complete  description  effectiveness  of 
this  iterative  algorithm. 

We  have  noted  before  that  the  linear  system  obtained  from  the  meshless 
method  is  consistent  even  if  the  stiffness  matrix  is  positive  semi-definite.  In 
this  situation,  a  solver  based  on  conjugate  gradient  method  can  also  be  used. 
The  convergence  in  this  situation  is  similar  to  the  convergence  of  conjugate 
gradient  method  when  applied  to  solve  the  linear  system  obtained  from  the 
standard  finite  element  method.  The  multigrid  method  is  not  directly  applicable 
to  the  linear  system  when  the  stiffness  matrix  is  positive  semi-definite,  since  the 
eigenfunctions  corresponding  to  the  zero  eigenvalue  of  the  stiffness  matrix  is 
global  and  oscillatory.  The  same  is  also  true  when  the  the  stiffness  matrix  is 
“almost”  singular.  However,  a  special  version  of  multigrid  method  was  proposed 
in  [76]  and  [44],  when  the  underlying  partition  of  unity  was  reproducing  of  order 
fc  =  0. 

4.  A-posteriori  error  estimation,  adaptivity,  and  computation  of 
data  of  interest 

The  rigorous  theory  of  a-posteriori  error  estimation  was  developed  in  [17]  and 
other  estimates,  based  on  various  averaging,  were  also  used.  These  estimators 
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can  be  used  as  error  indicators  for  adaptivity  purposes.  For  adaptive  approaches 
in  meshless  method,  we  refer  the  reader  to  [76]  and  [23]. 

We  finally  mention  that  programming  the  meshless  method  is  an  important 
issue  and  it  requires  specific  concepts.  For  this  aspect  of  the  meshless  method, 
we  refer  to  [32,  45,  82]  and  [83]. 

9  Examples 

Meshless  methods  have  been  applied  to  linear  and  non  linear  elliptic  problems, 
as  well  as  to  problems  related  to  other  differential  equations;  we  refer  to  [56]. 
However,  it  is  essential  to  characterize  the  types  of  problems  where  this  method 
is,  or  could  be,  superior  to  standard  methods  ([21]). 

In  this  paper,  we  address  only  the  application  of  the  meshless  method  on 
a  class  of  linear,  elliptic  problems.  As  stated  in  the  introduction,  one  of  the 
main  advantages  of  meshless  methods  is  that  it  avoids  mesh  generation.  This 
is  essential  when  the  domain  in  complex.  Another  advantage  of  this  method  is 
that  it  allows  the  use  of  various  “special”  local  shape  functions  to  improve  the 
accuracy. 

The  Generalized  Finite  Element  Method  (GFEM)  was  elaborately  discussed 
in  [83]  and  it  was  shown  that  the  method  is  effective.  Three  types  of  meshes 
with  successive  refinements  were  used  in  that  paper  and  we  present  one  of  these 
meshes  in  Figure  9.1.  This  is  a  simple  finite  element  mesh  and  it  does  not  reflect 
the  geometry  of  the  underlying  domain.  Then,  using  the  linear  finite  elements 
as  partition  of  unity  and  special  functions  for  local  approximation,  an  improve¬ 
ment  in  the  rate  of  convergence  was  achieved.  Detailed  numerical  data,  with 
comments  on  various  aspects  of  the  method,  e.g.,  numerical  integration,  etc. 
were  presented  in  [83].  We  note  however,  that  though  the  domain  considered 
in  this  example  {i.e.,  the  domain  in  Figure  9.1)  was  simple,  and  classical  fi¬ 
nite  element  method  with  mesh  refinement  or  an  adaptive  procedure  could  have 
been  used,  the  analysis  and  data  presented  in  [83]  clearly  shows  the  scope  and 
potential  of  GFEM. 


Figure  9.1:  A  mesh  used  in  [83]  for  the  constructing  of  partition  of  unity 
in  the  context  of  GFEM  to  solve  a  problem  posed  on  the  domain. 
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As  mentioned  above,  the  power  of  GFEM  lies  in  handling  problems  where 
the  underlying  domain  has  complex  geometry.  Three  types  complex  domains, 
shown  in  Figure  9.2,  were  analyzed  by  the  GFEM  in  [83].  Another  complex 
domain  with  fibers,  analyzed  in  the  same  paper,  is  shown  in  Figure  9.3  where 
the  fiber  distribution  was  taken  from  [10].  To  construct  finite  element  meshes 
for  these  domains  is  very  complex  and  nearly  impossible.  “Handbook”  problems 
that  characterize  the  local  behavior  of  the  approximated  solution  {e.g.,  in  the 
neighborhood  of  a  crack,  fibers,  etc.)  were  used  to  construct  special  shape 
functions  for  these  problems. 


Figure  9.2:  The  three  types  of  domains  analyzed  by  GFEM  in  [83]. 

The  GFEM  has  an  advantage  in  dealing  problems  with  singularities  (in  the 
neighborhood  of  geometric  edges)  in  3D.  When  the  basic  finite  element  tetra¬ 
hedral  mesh  is  used  in  such  problems,  it  is  well  known  that  classical  edge  re¬ 
finement  is  cumbersome.  This  problem  was  handled  using  GFEM  in  [36] ,  where 
a  refinement  by  special  functions,  at  positions  indicated  by  an  error  indica¬ 
tor,  was  performed.  GFEM  was  also  used  to  handle  difficulties  stemming  from 
orthotropic  problems  in  [35]. 
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Figure  9.3:  The  problem  of  fiber  composite  type  analyzed  in  [83]. 


There  are  other  types  of  problems  where  the  GFEM  is  quite  effective.  They 
include  multisite  problems,  where  many  crack  configuration  are  present,  and 
crack  propagation  problems,  where  the  geometry  of  the  domain  changes.  The 
GFEM  could  be  used  in  such  problems  by  considering  local  approximating 
spaces  consisting  of  functions  that  are  discontinuous  over  the  cracks  and  in¬ 
cluding  a  singular  function  with  respect  to  the  tip  of  the  crack  into  the  local 
approximating  space.  Then  the  propagation  of  the  crack  is  computed  (using 
stress  intensity  factors);  the  old  singular  function  in  the  approximating  space 
is  replaced  by  a  new  singular  function,  new  discontinuous  functions  are  added 
in  the  same  space,  and  the  process  of  computing  the  propagation  of  crack  and 
changing  the  local  approximation  spaces  is  repeated.  In  this  process,  the  ma¬ 
trix  of  the  underlying  linear  system  at  a  particular  step  can  be  obtained  by 
augmenting  the  matrix  corresponding  to  the  previous  step  with  new  rows  and 
columns,  and  it  is  possible  to  solve  the  new  linear  system  using  Schur  comple¬ 
ment,  which  uses  the  previously  computed  data.  This  general  idea  was  used  in 
[56,  64]  without  using  the  previously  computed  data. 

The  GFEM  is  an  important  tool  in  approximating  solutions  of  elliptic  prob¬ 
lems  with  rough  coefficients  as  well  as  homogenization  problems.  We  mention 
that  the  usual  finite  element  method  may  give  extremely  poor  result  when  ap¬ 
plied  to  elliptic  problems  with  rough  coefficients,  as  shown  in  [18].  It  was  shown 
in  [16],  using  a  detailed  analysis,  that  GFEM  leads  to  the  same  rate  of  conver¬ 
gence  for  problems  with  rough  coefficients  as  when  the  coefficients  are  smooth. 

We  emphasize  that  in  this  paper,  we  considered  only  a  small  (but  important) 
family  of  problems.  We  showed  that  the  use  of  meshless  methods,  particularly 
the  use  of  GFEM,  on  these  problems  is  advantageous  in  comparison  to  the 
standard  finite  element  method.  Of  course,  there  are  other  types  of  problems, 
especially  non-linear  problems,  which  we  have  not  addressed  in  this  paper. 

Note:  The  authors  will  like  to  thank  Elsevier  Science  for  allowing  the  use 
of  Figures  9.1,  9.2,  and  9.3  in  this  paper.  They  were  published  in  the  journal 
Computer  Methods  in  Applied  Mechanics  and  Engineering. 

10  Some  Comments  and  Future  Challenges 

In  the  previous  sections  we  addressed  linear  elliptic  equations.  The  approxi¬ 
mation  theory  we  developed  is  obviously  usable  in  practically  all  variationally 
formulated  problem,  provided  the  inf-sup  condition  (the  BB  condition)  is  sat¬ 
isfied.  Meshless  methods  and  the  GFEM  can  be  directly  used  for  higher  order 
equations  because  there  are  no  difficulties  in  constructing  shape  functions  of  any 
regularity.  Adaptive  procedures,  shock  capturing,  etc.,  can  be  combined  with 
adaptive  construction  of  shape  functions.  A  deep  theoretical  understanding  of 
these  issues  is  still  in  the  future. 

For  non-coercive  problems,  the  question  of  the  inf-sup  condition  is  more 
delicate  and  plays  an  important  role.  Some  non-coercive  problems  have  been 
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successfully  treated.  We  mention,  for  example,  equations  whose  solutions  are 
oscillatory,  as  with  the  Helmholtz  equation.  Special  shape  functions  were  used  to 
capture  the  oscillatory  character  of  the  solution;  see,  e.g.  [52],  [62].  Although 
there  are  some  problems  with  round-off  error,  it  can  be  expected  that  these 
difficulties  can  be  overcome. 

Meshless  Methods  can  be  applied  to  nonlinear  as  well  as  linear  problems. 
The  Handbook  approach  for  constructing  shape  functions  is  important  for  both 
types  of  problems. 

The  Finite  Element  Method  is  presently  used  for  solving  a  wide  variety  of 
problems.  Meshless  methods,  especially  the  GFEM,  offer  many  new  opportu¬ 
nities.  Since  the  FEM  is  a  special  case  of  the  GFEM,  the  study  of  the  GFEM 
provides  the  potential  to  develop  enhanced  and  improved  methods.  Under  the 
umbrella  of  GFEM,  there  is  actually  a  family  of  particular  methods — depending 
on  the  approximating  functions  employed.  The  goal  in  creating  new  particular 
GFEM  is  to  obtain  methods  that  are  especially  effective  on  certain  classes  of 
problems.  To  address  this  issue  meaningfully  and  effectively,  further  mathemat¬ 
ical  study  of  meshless  methods  is  essential. 
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